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I.  INTRODUCTION 


A.  BACKGROUND 

During  the  Persian  Gulf  War  UAVs  (Unmanned  Aerial  Vehicles)  proved  useful  by 
providing  battlefield  surveillance,  targeting  information,  and  naval  gunfire  spotting. 

They  were  able  to  provide  the  battlefield  commander  with  up  to  the  minute  intelligence 
and  offered  the  battlefield  commander  an  over  the  horizon  intelligence  gathering  asset 
that  possessed  the  flexibility  that  he  required  to  respond  to  the  constantly  changing 
circumstances  and  priorities  of  the  battlefield.  UAVs  also  provided  a  medium  range 
aerial  reconnaissance  and  intelligence  gathering  capability  that  was  not  hindered  by  the 
long  lead  time  required  for  other  reconnaissance  assets.  Current  UAV  capabilities  are 
primarily  limited  to  optical  and  infra-red  surveillance,  but  the  UAV  would  be  very  useful 
in  electronic  intelligence  (ELINT)  gathering.  UAVs  could  remain  on  station  listening  to 
the  enemy's  electronic  emissions  for  long  periods  of  time,  deep  in  enemy  territory, 
without  risking  aircraft  and  pilots. 

One  application  of  electronic  intelligence  gathering  is  locating  radar  emitters  from 
time  difference  of  arrival  (TDOA)  observations.  A  time  domain  technique  is  presented 
here  that  uses  the  Kalman  filter  to  estimate  the  location  of  a  radar  emitter  from  the  time 
difference  of  arrival  for  a  burst  of  pulses  between  two  or  more  receivers,  and  time 
difference  of  arrival  of  individual  pulses  between  two  or  more  sensors.  The  background 
information  and  models  used  to  develop  the  Kalman  filter  are  presented,  and  the  resulting 
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Kalman  filter  algorithms  are  tested  for  a  variety  of  problem  configurations.  Conclusions 
are  drawn  about  the  usefulness  of  the  algorithms  and  recommendations  are  made  for 
further  testing  and  evaluation. 

B.  RADAR  EMITTER  SIGNAL  CHARACTERISTICS 

The  signal  characteristics  most  important  in  the  emitter  location  problem  are  the 
Pulse  Repetition  Interval  (PRI),  and  the  pattern  and  rate  at  which  the  emitter  scans  a 
target.  These  characteristics  determine  which  estimation  technique  will  be  most 
effective,  and  they  affect  the  ability  of  these  algorithms  to  accurately  estimate  the  location 
of  the  emitter. 

The  PRI  is  the  period  of  time  between  successive  pulses.  This  value  can  range  from 
tens  of  milliseconds  for  long  range  search  radars  to  microseconds  for  pulse  Doppler 
radars.  Range  in  many  radar  systems  is  determined  by  measuring  the  time  required  for  a 
pulse  to  strike  a  target  and  return,  and  the  unambiguous  range  of  a  radar  system  is 
determined  by  the  length  of  the  PRI.  The  longer  the  PRI  the  longer  the  unambiguous 
range  of  the  radar.  Similarly,  the  PRI  affects  the  allowable  separation  of  the  sensors  used 
to  detect  the  pulse  TDOA.  The  longer  the  PRI  the  greater  the  allowable  separation  of  the 
sensors. 

The  vast  majority  of  radar  systems  are  designed  to  scan  a  volume  of  space  to  look  for 
targets.  All  of  these  systems  use  a  search  mode  that  is  used  to  scan  the  volume  of  space 
around  the  emitter  in  a  regular  pattern.  Most  radar  emitters  spend  the  majority  of  their 
operation  in  this  mode  and  it  is  in  this  mode  that  most  emitters  can  be  passively  detected. 
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The  method  that  the  emitter  uses  to  scan  and  search  for  a  target  is  imp>ortant  in  the 
problem  of  estimating  its  location  from  its  emitted  signals.  Radar  systems  typically  scan 
their  main  beam  mechanically  or  electronically. 

Mechanically  scarmed  radar  systems  designed  for  target  location  typically  scan  in 
azimuth  in  a  circular  pattern.  Their  antennas  are  rotated  by  motors  and  actuators  in  a 
circular  manner,  and  are  designed  to  provide  360  degree  search  coverage.  Special  radar 
systems  designed  for  altitude  finding,  or  tracking  and  fire  control  will  not  necessarily 
scan  in  azimuth  alone.  Most  of  these  systems  use  electronically  scanned  antennas  to  steer 
the  main  beam  of  the  radar,  and  can  scan  in  elevation  as  well  as  azimuth.  These  systems 
can  have  very  complex  scan  patterns  that  can  be  changed  as  required.  This  thesis  will  be 
primarily  concerned  with  mechanically  scanned  radar  emitters  that  scan  circularly  in 
azimuth.  The  majority  of  the  emitters  encountered  will  be  of  this  type.  The  principles 
developed  to  locate  these  emitters  can  also  be  applied  to  emitters  with  more  complex  scan 
patterns. 

As  an  emitter  scans  its  main  beam  past  a  target,  the  amplitude  of  the  pulses  that 
strike  the  target  will  generally  vary  from  zero  to  a  maximum  that  occurs  near  the  center  of 
the  main  beam.  Figure  1  shows  the  typical  burst  of  pulses  detected  as  the  emitter  scans 
past  a  receiver.  The  frequency  of  these  bursts  is  a  function  of  the  scan  rate  of  the  emitter, 
and  the  width  of  the  burst  is  a  ftmction  of  the  scan  rate  of  the  emitter  and  the  beamwidth. 
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Figure  1.  Burst  of  pulses  from  a  scanning  emitter 

C.  EMITTER  LOCATION  PROBLEM  DESCRIPTIONS 

Two  basic  approaches  are  used  to  estimate  the  location  of  the  emitter  from  the 
TDOA  of  its  signal.  In  the  first  approach  the  TDOA  for  the  burst  of  pulses  detected  as 
the  emitter  scans  past  receivers  is  used  to  estimate  the  location  of  the  emitter.  The  second 
approach  uses  the  TDOA  for  individual  pulses  between  two  sensors  on  widely  separated 
platforms  to  estimate  the  location  of  the  emitter. 

1.  Burst  TDOA  Problem  Description 

The  first  TDOA  filtering  problem  assumes  that  multiple  receivers  are  used  and 
that  they  are  widely  separated.  The  TDOA  for  the  burst  of  pulses  detected  at  each 
receiver  are  measiued  and  used  to  estimate  the  location  of  the  emitter.  In  this  burst 

TDOA  filtering  problem,  the  following  assumptions  are  made; 

1 .  Enough  information  is  obtained  fix>m  the  received  signals  to  associate  TDOA 
observations  with  the  proper  emitter,  and  thus  the  multi-emitter  problem  is 
reduced  to  a  single  emitter  problem. 


2  The  receivers  are  referenced  to  a  common  ume  base  All  time  of  amvai 
measurements  are  referenced  to  this  common  time  base 

3  GlohaJ  Positioning  System  (GPS)  is  used  to  determine  the  locations  of  the 
receivers.  These  estimates  of  the  receiver  locations  are  assumed  to  be  exact  and 
errm  introduced  by  the  GPS  is  ignored. 

4.  The  emitter  scans  in  azimuth  at  a  constant  rate. 

5.  The  PRl  is  constant  for  the  entire  length  of  the  burst. 

The  receivers  may  be  mounted  on  UAVs,  aircraft,  ground  sites,  or  space  platforms.  For 
the  UAV  mounted  receivers,  the  desired  maximum  detection  range  is  assumed  to  be 
500  kilometers.  The  algorithms  developed  are  tested  for  ranges  of  500  kilometers  as  well 
as  shorter  ranges  of  1 50  kilometers  and  30  kilometers. 

2.  Single  Pulse  TDOA  Problem  Description 

The  pulse  TDOA  problem  estimates  the  location  of  the  emitter  by  measuring  the 
TDOA  for  an  individual  pulse  between  two  sensors.  The  sensors  are  spaced  close  enough 
so  that  the  TDOA  can  be  measured  for  a  single  pulse  without  ambiguity.  The  following 

assumptions  are  made  in  the  problem  development: 

1 .  Enough  information  can  be  obtained  from  the  received  signals  to  associate  TDOA 
observations  with  the  proper  emitter,  and  thus  the  multi-emitter  problem  is 
reduced  to  a  single  emitter  problem. 

2.  The  receiver  platforms  are  referenced  to  a  common  time  base.  All  time  of  arrival 
measurements  are  referenced  to  this  common  time  base. 

3.  Each  receiver  platform  has  two  sensors  and  the  TDOA  for  individual  pulses 
between  the  two  sensors  can  be  measured  accurately. 

4.  GPS  is  used  to  determine  the  locations  of  the  receiver  platforms.  These  estimates 
of  the  reiver  locations  are  assumed  to  be  exact  and  error  introduced  by  the  GPS 
is  ignored. 

The  algorithms  are  tested  for  ranges  of  500  kilometers  as  well  as  shorter  ranges  of  1 50 
kilometers  and  30  kilometers. 
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II.  TDOA  ESTIMATION 


A.  TIME  OF  ARRIVAL  ESTIMATION  AND  ERROR  MODELS 

The  Kalman  filter  can  be  used  to  estimate  the  location  of  the  emitter  from  TDOA 
observations.  Kalman  filter  theory,  however,  assumes  that  the  noise  corrupting 
observations  is  Gaussian  distributed.  Before  applying  the  Kalman  filter  to  the  emi  tier 
location  problem,  an  analysis  of  the  error  statistics  of  the  TDOA  observations  is 
necessary.  This  determines  if  the  noise  present  in  TDOA  observations  can  be 
approximated  as  Gaussian.  One  method  of  estimating  the  TDOA  for  emitter  bursts  and 
pulses  is  presented.  The  effect  of  noise  on  this  estimation  method's  accuracy  is  explored, 
and  the  probability  densities  calculated  are  compared  to  Gaussian  distributions. 

1.  Time  of  Arrival  Measurement 

TDOA  measurements  between  widely  separated  receivers  and  sensors  are  used  to 
estimate  the  location  of  the  emitter.  Since  the  receivers  need  to  be  widely  separated  to 
improve  the  accuracy  of  the  location  estimates  and  the  orthogonality  of  the  TDOA 
observations,  physically  linking  the  receivers  together  is  not  feasible.  Therefore,  it  is 
important  that  the  receivers  are  accurately  referenced  to  some  external  time  base.  An 
onboard  time  base  can  be  used  to  estimate  the  time  of  arrival  (TOA)  of  signals,  but  drift 
would  introduce  error  into  the  estimates  of  the  TOA  and  subsequently  the  TDOA 
observations.  To  minimize  drift,  internal  clocks  aboard  the  receivers  can  be  periodically 
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corrected  with  references  to  a  more  accurate  external  time  base.  If  the  receivers  are 
mounted  on  UAVs  equipped  with  GPS,  the  GPS  system  time  could  be  used  as  the 
external  time  base.  This  approach  is  assumed  for  the  balance  of  the  development. 

The  TDOA  observations  in  the  burst  TDOA  filtering  problem  are  on  the  order  of 
milliseconds  and  tens  of  milliseconds.  The  error  typically  present  in  the  GPS  system 
time  is  much  smaller  than  the  typical  TDOA  observations  and  synchronizing  all  of  the 
receivers  to  the  GPS  system  time  would  not  introduce  appreciable  error  into  the 
observations.  For  these  observations,  the  receivers  are  assumed  individually 
synchronized  to  the  GPS  system  clock. 

For  the  pulse  TDOA  observations,  the  length  of  the  typical  observations  are  a 
few  microseconds  or  less.  Using  the  GPS  system  time  to  estimate  the  time  of  arrival  of  a 
pulse  at  a  sensor  would  introduce  unacceptable  error  in  the  TDOA  observations.  This  is 
because  the  error  present  in  the  GPS  system  clock  would  be  on  the  same  order  as  the 
observation  it  is  being  used  to  measure.  To  avoid  ambiguity  in  the  pulse  TDOA 
observations,  the  sensors  must  be  placed  relatively  close  together.  For  the  balance  of  the 
developments  here  the  sensors  will  be  assumed  physically  linked  to  the  receiver  platform, 
and  the  TOA  observations  for  both  sensors  will  be  referenced  to  the  same  on  board  time 
base. 

2.  Pulse  TOA  Estimation 

Envelope  detection  of  the  pulses  received  by  the  sensors  is  assumed.  The 
envelopes  of  the  received  pulses  are  sampled  and  digitally  encoded.  An  estimate  of  the 
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TOA  for  a  pulse  can  be  obtained  from  the  time  at  which  a  sample  is  detected  above  a  set 
threshold.  This  method  of  estimating  the  TOA  is  highly  sensitive  to  the  sampling  rate 
and  noise  present  in  the  sampled  pulse.  Spurious  noise  that  is  above  threshold  gives  false 
TOA  indications.  Algorithms  can  be  developed  to  test  for  false  indications  by  examining 
multiple  samples  and  announcing  the  arrival  of  a  pulse  only  if  a  specified  number  of 
above  threshold  samples  are  detected.  A  more  accurate  estimate  of  the  TOA  can  be 
obtained  by  integrating  the  samples  of  the  pulse  and  calculating  the  TOA  of  the  centroid 
of  the  pulse.  The  improvements  gained  from  integration  are  similar  to  the  improvements 
obtained  in  radar  systems  when  pulse  integration  is  used.  All  of  the  information  present 
in  the  pulse  is  utilized  and  the  effective  signal  to  noise  ratio  is  increased. 

The  time  of  arrival  of  the  centroid  of  a  sampled  pulse  can  be  found  from 
equation  (1): 


TOA  = 


£tOc)z(k) 

ks\ 


(1) 


Where  TOA  =  the  Time  Of  Arrival  (TOA)  of  the  pulse  centroid, 

z(k)  =  the  magnitude  of  the  kth  sample, 
t(k)  =  the  time  at  which  the  kth  sample  was  taken,  and 
N  =  the  number  of  samples  in  the  pulse. 

In  the  following  development  of  the  error  present  in  the  estimate  of  the  pulse 


TOA,  detected  pulses  are  assumed  to  have  the  shape  shown  Figure  2.  To  these  pulses. 


white,  Gaussian  distributed  noise  is  added  and  its  effect  on  the  TOA  is  calculated. 
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Figure  2.  Pulse  envelope  and  sampled  pulse 


The  sampling  interval  is  assumed  to  be  much  smaller  than  the  pulse  width  so  a 
sufficient  number  of  samples  of  the  pulse  are  taken  and  the  sampled  pulse  is  a  reasonable 
representation  of  the  original  pulse.  The  sampling  is  assumed  to  be  triggered  by,  and 
aligned  with,  the  rising  edge  of  the  pulse.  This  assumption  ignores  the  error  introduced  in 
the  TOA  when  the  first  samples  of  the  pulse  are  taken  at  some  random  point  within  the 
pulse. 
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Each  sample  of  the  pulse  is  considered  a  random  variable  as  shown  in 
equation  (2). 

2(k)  =  Zo(k)  +  v(k)  (2) 

Where  z(k)  =  the  random  variable  representing  the  amplitude  of  the  kth  sample, 

Zo(k)  =  the  deterministic  variable  that  is  the  true  amplitude  of  the  kth 
sample,  and 

v(k)  =  the  white  Gaussian  random  variable  that  represents  the  noise 
present  in  the  amplitude  of  the  kth  sample.  The  noise  has  zero 
mean  and  a  variance  of  V,  i.e..  E[v(k)^]  *  V. 

Since  each  sample  of  the  pulse  is  modeled  as  a  random  variable,  the  time  of  arrival  of  the 
centroid  of  the  pulse,  calculated  in  equation  (1),  is  also  a  random  variable.  A  detailed 
derivation  of  the  probability  density  of  the  pulse  TOA  is  presented  in  Appendix  A. 

The  peak  signal  to  noise  ratio  is  used  as  a  means  of  representing  the  noise  power 
present  in  the  samples  of  the  pulse.  As  shows  in  Figure  2,  the  peak  amplitude  of  the 
pulse  is  assumed  to  be  one.  The  peak  signal  to  noise  power  ratio  is  calculated  from  the 
following  equations: 

f  Speik"!  _  E[Zo(k)^] 

^  N  J  E[v(k)2]  ^ 

Where  z^Ck)  -  the  true  amplitude  of  the  peak  sample  of  the  pulse,  and 

v(k)  =  the  zero  mean  white  noise  sample  with  variance  V,  and 
E[  ]  =  the  expectation  operator. 

In  decibels  the  peak  signal  to  noise  power  ratio  is  given  as: 
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Since  the  peak  amplitude  of  the  pulse  is  one  and  the  noise  v(k)  is  a  random  variable  with 
zero  mean  and  a  variance  V,  the  peak  signal  to  noise  power  ratio  reduces  to: 

,5, 

The  MATLAB  program  Puldist.m,  listed  in  Appendix  E,  calculates  the 
probability  densitiy  functions  plotted  in  Figure  3.  These  curves  depict  the  probability  that 
the  TOA  of  a  pulse,  calculated  using  equation  (1),  will  fall  within  a  specified  range.  The 
mean  and  standard  deviation  of  the  probability  density  fimctions  are  also  calculated  and 
Gaussian  distributed  density  fimctions  having  the  same  mean  and  standard  deviation  are 
superimposed  on  the  plots. 


Figure  3.  Pulse  TOA  probability  densities  vs.  Gaussian 
approximations,  pulsewidth  1.0  psec 
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The  Gaussian  density  functions  are  nearly  indistinguishable  from  the  actual 
probability  density  functions,  so  the  pulse  TOA  can  be  considered  Gaussian  distributed  in 
the  Kalman  filter  developments.  In  Appendix  A,  a  complete  analysis  is  conducted  of  the 
effects  of  peak  signal  to  noise  ratio  and  sampling  rate  on  the  probability  density  function 
of  the  pulse  TOA. 

3.  Burst  TOA  Estimation 

As  the  main  beam  of  an  emitter  sweeps  past  a  receiver,  the  receiver  detects  a 
burst  of  pulses.  The  number  of  pulses  received  is  a  function  of  the  beam  width  and  scan 
rate  of  the  emitter.  The  amplitude  of  the  pulses  change  as  the  main  beam  of  the  emitter 
sweeps  past  the  receiver.  A  burst  typical  for  circularly  scanning  emitters  is  shown  in 
Figure  1 .  The  envelope  of  the  burst  can  be  approximated  with  the  upper  half  of  a  sine 
wave. 

x(t)  =  sm[^®t]  (6) 

Where  x(t)  =  the  amplitude  of  the  envelope  of  the  burst, 

BW  =  the  beamwidth  of  the  main  beam  of  the  emitter  in  radians, 

SR  =  the  scan  rate  of  the  emitter  in  radians/sec,  and 
t  =  the  time  in  seconds. 

An  estimate  of  the  burst  TOA  can  be  obtained  from  the  time  of  arrival  for  the 
first  pulse  of  a  burst,  but  this  method  of  estimating  the  TOA  assumes  that  for  each 
receiver  the  same  pulse  in  the  burst  will  be  detected  first.  Because  of  lower  signal 
strength,  more  distant  receivers  may  not  detect  some  pulses  as  the  start  of  the  burst  and 
the  estimate  of  the  TOA  will  be  in  error  by  some  multiple  of  the  PRI.  Estimating  the 
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centroid  of  the  envelope  of  the  burst  and  its  time  of  arrival  yields  a  better  estimate  of  the 
TOA  because  the  shape  of  the  envelope  detected  by  all  of  the  receivers  is  the  same,  and 
only  the  amplitude  of  the  envelope  changes  as  a  function  of  distance  from  the  emitter. 

An  estimate  of  the  TOA  of  the  burst  envelope  is  found  by  integrating  the  sampled 
pulses  and  their  time  of  arrival  over  the  length  of  the  burst.  If  enough  pulses  are  present 
in  the  burst,  the  centroid  of  the  pulses  lies  sufficiently  close  to  the  centroid  of  the 
envelope.  Integration  of  the  pulses  also  reduces  the  effect  of  the  noise  by  averaging  its 
effect  over  a  longer  period  of  time,  and  effectively  improving  the  signal  to  noise  ratio  by 
ignoring  the  interpulse  periods  where  only  noise  is  present. 

Since  the  burst  is  considered  a  collection  of  sampled  pulses,  equation  (1 )  is  used 
to  estimate  the  time  of  arrival  of  the  centroid  of  the  burst.  The  equations  and  algorithms 
used  to  calculate  the  probability  density  function,  mean,  and  standard  deviation  of  the 
pulse  TOA  are  used  to  calculate  these  properties  for  the  burst  TOA. 

The  MATLAB  program  Burdist.m,  listed  in  appendix  E,  generates  the  burst  in 
Figure  4  and  calculates  its  probability  density  for  peak  signal  to  noise  ratios  of  IS  and 
25  dB. 

The  burst  has  the  following  characteristics: 

Burst  Length:  2.0  milliseconds 

Pulsewidth  1 .00  microseconds. 

Pulse  Repetition  Frequency  6000  Hz, 

Maximum  Amplitude  1 .00 
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Figure  4.  A  burst  of  pulses  with  a  6.0  kHz  PRF  and 
a  1.0  microsecond  pulsewidth 


The  envelope  of  the  burst  is  calculated  using  equation  (6),  and  the  peak  signal  to 
noise  ratios  are  calculated  using  equations  (4)  and  (5).  The  burst  in  Figure  4  assumes  that 
the  pulses  are  spaced  evenly  throughout  the  envelope  of  the  burst.  If  this  is  the  case,  the 
centroid  of  the  pulses  coincide  with  the  true  centroid  of  the  envelope.  For  emitters  with 
long  PRIs  where  few  pulses  are  present  in  the  burst,  this  assumption  may  not  be  valid, 
and  the  error  introduced  could  be  significant.  In  all  further  developments  emitter 
characteristics  are  chosen  so  that  this  error  is  insignificant. 

For  the  probability  density  function  of  the  burst  shown  in  Figure  4,  the  mean  and 
standard  deviation  are  calculated,  and  Gaussian  distributions  having  the  same  mean  and 
standard  deviation  are  simultaneously  plotted.  Figure  5  is  a  plot  of  the  burst  TOA 
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probability  density  function  and  the  Gaussian  approximation  for  a  peak  signal  to  noise 
ratio  of  1 5  dB.  Figure  6  is  the  plot  of  the  burst  TOA  probability  density  function  and 
Gaussian  approximation  for  a  peak  signal  to  noise  ratio  of  25  dB. 


Figure  5.  Burst  TOA  probability  density  vs.  Gaussian  approximation, 
2.0  millisecond  burst,  6.0  kHz  PRF  and  1.0  microsecond  pulsewidth 


15 


U| _ I _ y  _ i, _ I _ I 

0.9  0.95  1.0  1.05  1.1 

Centroid  Time  of  Arrival  (miUiseconcIs) 


Figure  6.  Burst  TOA  probability  density’  vs.  Gaussian  approximation, 

2.0  millisecond  burst,  6,0  kHz  PRF,  and  1.0  microsecond  pulsewidth 

The  burst  TOA  probability  density  functions  are  nearly  indistinguishable  from 
the  Gaussian  probability  density'  function.  The  burst  TOA  can  be  considered  Gaussian 
distributed  in  the  Kalman  filter  development.  In  Appendix  A  a  more  complete  analysis  is 
conducted  of  the  effects  of  peak  signal  to  noise  ratio,  pulse  width,  and  PRF  on  the 
probability  density  functions  of  the  burst  TOA. 


B.  BURST  AND  PULSE  TDOA 

As  demonstrated  in  the  previous  sections,  the  TOA  of  pulse  or  burst  can  be 
approximated  as  a  Gaussian  random  variable.  If  the  TOA  for  each  burst  or  pulse  is 
considered  independent,  the  TDOA  can  be  calculated  using  the  following  equation; 
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TDOA  =  Za  -  Zb 


(7) 


Where  2^=  TOA  of  burst  or  pulse  A,  and 
Zb  =  TOA  ofburst  or  pulse  B. 

Time  of  arrival  A  and  B  can  be  bursts  generated  by  an  emitter  scanning  two  different 
receivers,  or  the  bursts  generated  by  different  scans  of  the  emitter  past  the  same  receiver. 
Time  of  arrival  A  and  B  can  also  be  the  time  of  arrival  for  two  different  pulses  at  the 
same  sensor  or  the  time  of  arrival  of  a  single  pulse  at  two  different  sensors. 

If  the  TOA  is  expressed  as  the  sum  of  a  constant  value  uliich  corresponds  to  the  true 
TOA,  and  a  zero  mean  Gaussian  random  variable,  then  the  TOA  is: 


Za=ZoA+Va 


(8) 


Where  Za=  the  observed  time  of  arrival, 

ZoA=  the  true  time  of  arrival,  and 

Va  =  the  noise  present  in  the  observed  time  of  arrival. 

This  noise  has  zero  mean  and  variance  Va. 

The  noise  term  Va  is  considered  to  be  white  noise,  with  zero  mean  and  uncorrelated 
increments.  Therefore,  the  expected  value,  E[VaVb],  is  zero.  The  mean  of  the  estimate  of 
the  TDOA  are: 


Ptdoa  =  E[za-zb  ]  =  E[Zoa-Zob  ]  +  E[va  ]-E[vb  ]  (9) 


PTDOA  =  ZoA  -  ZoB 


(10) 


Therefore,  the  expected  value  of  the  TDOA  is  the  difference  between  the  true  times  of 
arrival  for  the  two  bursts,  i.e.  the  true  TDOA. 

The  variance  of  the  estimate  of  the  TDOA  is: 
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^TOOA  ~  E[  (TDOA  -  ^tdoa)^  ] 


(11) 


<y™A=E[(VA-VB)M  (12) 

Since  the  increments  of  the  noise  are  independent,  and  E[vaVb]  =  0,  the  resulting  variance 
of  the  TDOA  estimate  is: 

o™A  =  E[vi  +  v|]  =  V*  +  VB  (13) 

Therefore,  the  TDOA  of  a  burst  or  a  pulse  is  a  Gaussian  random  variable  with  a 
mean  value  equal  to  the  true  TDOA,  and  a  variance  equal  to  the  sum  of  the  variances  of 
the  burst  or  pulse  time  of  arrival  estimates. 
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ni.  EMITTER  LOCATION 


The  TDOA  for  a  radar  signal  is  the  difference  in  the  amount  of  time  that  a  radar 
signal  takes  to  reach  two  receivers.  This  research  utilizes  the  time  difference  of  arrival 
for  bursts  of  pulses  between  two  widely  separated  receivers  and  the  TDOA  for  individual 
pulses  between  two  closely  spaced  sensors  to  estimate  the  location  of  an  emitter.  Both  of 
these  TDOAs  are  iunctions  of  the  emitter  and  receiver  locations  and  the  emitter 
characteristics.  In  the  approaches  pursued  here  only  the  two  dimensional  problem  is 
considered,  and  the  effect  of  altitude  is  ignored. 

A.  BURST  TIME  DIFFERENCE  OF  ARRIVAL  (TDOA) 

1.  Problem  Geometry 

The  TDOA  of  a  radar  signal  burst  for  two  widely  separated  receivers  is  primarily 
due  to  the  amount  of  time  required  for  the  emitter  to  scan  the  two  receivers.  If  the 
assumption  is  made  that  the  emitter  scans  only  in  azimuth,  then  the  TDOA  is  only  a 
function  of  the  scan  rate  of  the  emitter  and  the  angle  formed  between  the  emitter  and  the 
two  receivers.  The  scan  rate  is  the  angular  rate  of  rotation  of  the  emitter,  and  the  angle 
formed  by  the  emitter  and  the  two  receivers  can  be  determined  from  the  coordinates  of 
the  emitter  and  receivers.  If  the  scan  rate  of  the  emitter  is  assiuned  constant  and  known 
a  priori  or  from  some  other  estimation  method,  then  the  TDOA  becomes  a  function  of 
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only  the  location  of  the  emitter  and  the  receivers.  Figure  7  is  the  geometric  relationship 
between  the  receivers,  and  the  emitter. 


Figure?,  Burst  TOO  A  geometry 


Since  the  scan  rate  is  known,  the  TOOA  for  an  emitter  scanning  two  receivers  is 
only  function  of  the  angle  6  formed  by  the  emitters  and  the  receivers.  This  angle  can  be 
found  by  applying  vector  geometry.  If  the  positions  of  the  receivers  are  known,  the  cross 
product  of  the  position  vectors  from  the  two  receivers  to  the  target  can  be  used  to  express 
the  angle  formed  by  the  two  receivers  and  the  emitter  in  terms  of  the  unknown  quantities, 
the  X  and  y  coordinates  of  the  emitter.  The  angle  0  is  calculated  from  the  following 
equations: 
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RaXRb  =  iRaliRbI  sin(0) 


(14) 


M  A  A 

Ra  =  (xt  -  xa)i  +  (yt  -  ya)j 


(15) 


Rb  =  (xt  -  xb)i  +  (yt  -  yb)j 


(16) 


RaXRb  = 


A  A  A 

i  j  k 
(xt-xa)  (yt-ya)  0 
(xt-xb)  (yt-yb)  0 


(17) 


sin(0)  =  (xt-xa)(yt-yb)-(yt-yaXxt-xb) 

[(xt  -  xa)^  +  (yt  -  ya)^]  ‘^[(xt  -  xb)^  +  (yt  -  yb)^] 

Where  Ra  =  The  position  vector  from  receiver  A  to  the  emitter, 
Rb  =  the  position  vector  from  receiver  B  to  the  emitter, 

0  =  the  angle  formed  by  the  receivers  and  the  emitter, 
xt,  yt  =  the  X  and  y  coordinate  of  the  emitter, 

xa, ya  *  the  x  and  y  coordinate  of  receiver  A, 

xb, yb  *  the  X  and  y  coordinate  of  receiver  B,  and 
SR  ==  the  scan  rate  of  the  antenna  in  rad/sec. 

If  0  is  small,  the  TDOA  is  found  from  the  following  equation. 


(18) 


TDOA  = 

SR  SR 


(xt-xa)(yt-yb)-(yt-yaXxt-xb) 

L  [(xt  -  xa)^  +  (yt  -  ya)^J  ''^[(xt  -  xb)^  +  (yt  -  yb)^]  J 


(19) 


2.  Loci  of  Constant  TDOA 

From  "quation  (19)  the  burst  TDOA  for  two  widely  separated  receivers  being 
scanned  by  a  constantly  rotating  emitter  is  directly  proportional  to  the  angle  formed  by 
the  two  receivers  and  the  emitter.  In  Appendix  B  an  approximate  relationship  is 
developed  to  plot  the  loci  of  points  where  an  emitter  could  lie  and  produce  a  specific 
TDOA  observation.  These  loci  of  constant  TDOA  are  shown  in  Figure  8. 
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Figure  8.  Loci  of  constant  TDOA  for  an  emitter  scanning  at  a  constant  rate, 
scan  rate  360  degrees/sec,  receiver  separation  1.0  kilometers 

These  loci  demonstrate  that  for  a  single  TDOA  observation  there  exist  an  infinite 
number  of  possible  location  for  the  emitter.  Another  observation  or  measurement  is 
required  to  uniquely  estimate  the  location  of  the  emitter. 

If  more  than  two  receivers  are  used,  a  TDOA  observation  is  available  for  each 
combination  of  receiver  pairs.  Plotting  the  loci  of  constant  TDOA  for  each  of  the 
receiver  pairs,  multiple  intersections  occur,  but  all  of  the  loci  intersect  at  a  common  point, 
the  location  of  the  emitter.  This  is  shown  in  Figure  9. 
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Figure  9.  Loci  of  constant  TDOA  for  an  emitter  scanning  at  a  constant  rate, 
scan  rate  360  degrees/second,  receiver  locations  as  shown 

3.  Orthogonality  of  TDOA  Observations 

The  performance  of  the  Kalman  filter  varies  based  upon  the  quality  of  the 
observations  available.  Generally,  the  observations  with  the  highest  orthogonality 
perform  the  best.  If  data  processing  time  or  capability  is  limited,  it  is  advantageous  to  use 
the  observations  with  the  highest  orthogonality.  The  loci  of  constant  TDOA  are  one  way 
of  visualizing  and  calculating  the  orthogonality  of  the  observations. 

A  measure  of  the  orthogonality  of  the  TDOA  observations  is  obtained  from  the 
dot  product  of  the  unit  vectors  tangent  to  the  loci  of  constant  TDOA.  These  vectors  are 
calculated  from  the  slope  of  the  line  tangent  to  the  loci  at  the  estimated  location  of  the 
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emitter.  Once  the  unit  vectors  tangent  to  the  loci  are  found,  their  dot  product  yields  the 
cosine  of  the  angle  formed  by  the  two  loci.  The  equation  for  the  dot  product  is  given 
below: 


Rti»Rt2=cos0  (20) 

Where  Rt,  =  The  unit  vector  tangent  to  the  loci  munber  1, 

Rt}  =  the  unit  vector  tangent  to  loci  number  2,  and 
0  =  the  angle  formed  by  the  two  loci. 

The  cosine  of  the  angle  formed  by  the  loci  is  an  indication  of  the  orthogonality  of 
the  loci.  If  the  cosine  of  0  is  close  to  one  then  the  loci  are  nearly  parallel,  and  if  the 
cosine  of  0  is  zero  the  loci  are  orthogonal  and  lie  at  right  angles  to  each  other. 

The  unit  vectors  tangent  to  the  loci  are  found  using  the  following  equation: 

Rt  = -7=====-X  + -T^rrry 

JT+n?  Vl  +m^  ^ 

Where  Rt  =  the  unit  vector  tangent  to  the  loci,  and 

m  =  the  slope  of  the  loci,  dyt/dxt  at  the  emitter. 

The  equations  for  the  slope  of  the  loci  at  an  arbitrary  point  are  developed  in 


Appendix  C  and  presented  here: 


dyt  _  [(xt  -  xa)^sin  y  -  2(xt  -  xa)(yt  -  ya)cos  v|/  -  (yt  -  ya)^sin  y] 
dxt  [(yt  -  ya)^cos  vj/  -  2(xt  -  xa)(yt  -  ya)sin  v}/  -  (xt  -  xa)^cos  v] 

Where  xt,  yt  =  the  x  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  X  and  y  coordinates  of  receiver  A, 

xb,  yb  =  the  X  and  y  coordinates  of  receiver  B,  and 
vf;  =  the  bearing  from  receiver  A  to  receiver  B. 

The  bearing  from  receiver  A  to  receiver  B  is  calculated  from  the  following  equation: 


V  =  arctan 


fyb-ya' 

Lxb-xa. 


4.  Ambiguity  in  TDOA  Observations 


Ambiguity  in  the  burst  TDOA  observations  is  present  if  another  burst  is  received 
by  one  of  the  receivers  before  all  of  the  receivers  detect  the  first  burst.  In  the  burst 
TDOA  problem  each  of  the  receivers  is  scaimed  only  once  per  rotation  of  the  emitter. 
Ambiguity  is  possible  if  the  distance  separation  between  the  receivers  causes  a  burst  fi-om 
the  next  scan  of  the  emitter  to  be  received  by  a  closer  receiver  before  a  distant  receiver 
detects  the  burst  fi-om  the  first  scan.  This  occurs  if  the  time  required  for  the  signal  to 
travel  firom  one  receiver  to  the  next  is  greater  than  the  time  required  for  the  emitter  to 
complete  one  scan.  This  is  unlikely.  In  even  the  fastest  scanning  emitters  the  scan  rate  is 
unlikely  to  be  greater  than  10  revolutions  per  second.  For  an  emitter  with  a  scan  rate  of 
10  revolutions  per  second,  the  separation  required  to  introduce  ambiguity  in  the  burst 
TDOA  is  30,000  kilometers. 

B.  SINGLE  PULSE  TIME  DIFFERENCE  OF  ARRIVAL  (TDOA) 

1.  Problem  Geometry 

If  two  sensors  are  illuminated  at  the  same  time  by  the  beam  of  an  emitter  sending 
out  pulsed  signals,  the  pulse  TDOA  is  the  difference  in  the  amount  of  time  required  for  a 
single  pulse  to  reach  the  two  sensors.  Since  electromagnetic  waves  travel  at 
approximately  the  speed  of  light,  the  TDOA  is  a  fimction  of  the  difference  in  the  distance 
that  the  pulse  must  travel.  The  difference  in  distance  the  pulse  must  travel  is  calculated 
from  the  geometry  of  the  problem.  Figure  10  is  a  diagram  of  the  geometry  of  the 
problem. 


25 


Figure  10.  Single  pulse  TDOA  geometry 


The  time  difference  of  arrival  for  a  pulse  at  sensors  A  and  B  is  calculated  from 
equation  (24). 


[(xt  -  xa)^  +  (yt  -  ya)^  ]  -  [(xt  -  xb)^  +  (yt  -  yb)^  ] 

TDOA  = - * — ^ - - — 

Where  TDOA  =  the  time  difference  of  arrival, 

xt,  yt  =  the  X  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  x  and  y  coordinates  of  sensor  A, 

xb,  yb  =  the  X  and  y  coordinates  of  sensor  B, 
c  =  the  speed  of  light. 


(24) 
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2.  Loci  of  Constant  TDOA 


Given  a  pulse  TDOA  observation  and  sensor  geometry,  the  possible  location  of 
the  emitter  is  not  unique  and  a  locus  of  possible  emitter  locations  exists.  Figure  1 1  is  a 
plot  of  equation  (24)  for  constant  values  of  TDOA  and  the  sensor  configuration  shown. 
The  sensors  are  separated  by  500  meters.  For  the  TDOA  observations  shown,  the 
location  of  the  emitter  lies  anywhere  along  these  hyperbolic  loci. 


Figure  11.  Loci  of  constant  pulse  TDOA, 
sensor  separation  500  meters 


If  more  than  two  sensors  are  used  it  is  possible  to  determine  a  unique  estimate  of 
the  emitter's  location.  The  scenario  shown  in  Figure  12  uses  two  widely  separated 
receiver  platforms  with  two  sensors  attached  to  each  of  the  receiver  platforms.  With  this 
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configuration,  the  loci  of  constant  TDOA  intersect  at  a  single  point  and  a  unique  estimate 
of  the  emitter's  location  can  be  obtained  from  the  pulse  TDOA  observations. 
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Figure  12.  Loci  of  constant  TDOA  for  two  receiver  platforms, 
sensor  separation  500  meters 


3.  Orthogonality  of  the  Loci  of  Constant  TDOA 

As  previously  developed  in  the  burst  TDOA  problem  the  orthogonality  of  the 
TDOA  observations  can  be  related  to  the  orthogonality  of  the  loci  of  constant  TDOA. 
The  dot  product  of  the  unit  vectors  tangent  to  the  loci  give  the  cosine  of  the  angle 
between  the  loci  and  a  good  indication  of  the  orthogonality  of  the  TDOA  observations. 
The  unit  vectors  tangent  to  the  loci  are  calculated  using  equation  (21).  The  equations  for 
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the  slope  of  the  loci  of  constant  TDOA  are  developed  in  Appendix  C.  The  slope  of  the 
loci,  dxt/dyt,  is  given  by  the  following  equation: 


dyt  _  (yt  -  ya)^eos  ly  -  (xt  -  xaXyt  -  ya)sin  y 
dxt  (xt-xa)(yt-ya)cosv|/-(xt-xa)^sLn\|/ 

Where  xt,  yt  =  the  x  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  X  and  y  coordinates  of  sensor  A, 

xb,  yb  =  the  X  and  y  coordinates  of  sensor  B,  and 
v(/  =  the  bearing  from  sensor  A  to  sensor  B. 

The  bearing  frt>m  sensor  A  to  sensor  B  is  calculated  from  the  following: 


=  arctan 


(yb-ya) 

(xb-xa) 


(25) 


(26) 


4.  Ambiguity  in  Single  Pulse  TDOA  Observations 

Ambiguity  in  pulse  TDOA  observations  develops  if  a  second  pulse  is  detected  by 
one  of  the  sensors  before  both  of  the  sensors  detect  the  first  pulse.  This  can  occur  if  the 
amount  of  time  required  for  a  pulse  to  travel  from  one  sensor  to  the  next  is  longer  than  the 
PRI.  Therefore,  the  separation  of  sensors  determines  the  highest  PRF  that  signals  can 
have  and  still  be  detected  without  ambiguity. 

PRFm«  = 

PRI  =  Rab  (27) 

Where  Rab  =  the  distance  between  sensor  A  and  sensor  B,  and 

c  =  the  speed  of  light. 

For  sensors  separated  by  a  distance  of  500  meters  the  smallest  PRI  an  emitter  could  have 
and  be  detectable  without  ambiguity  is  1 .66  microseconds.  This  corresponds  to  a  PRF  of 
approximately  600  kHz. 
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Estimates  of  the  PRI  of  the  emitter  can  be  calculated  from  the  time  interval 


between  pulses  at  each  of  the  sensors.  Therefore,  the  PRI  can  be  estimated  independently 
and  used  to  test  the  TDOA  observations  for  ambiguity. 
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rv.  KALMAN  FILTERING 


A.  THE  EXTENDED  KALMAN  FILTER 

There  are  several  methods  of  linearizing  a  nonlinear  problem  to  allow  the  Kalman 
filter  to  be  used.  The  problem  itself  can  be  linearized  by  expanding  the  equations  of  state 
in  a  Taylor  series  around  some  norm.  This  approach  is  often  used  to  derive  equations  of 
state  in  aircraft  stability  and  control  problems.  The  nonlinear  equations  of  state  are 
expanded  and  derived  for  pertubations  around  some  equilibrium  trim  point.  As  long  as 
the  pertubations  from  the  trim  point  are  within  acceptable  limits  this  linear  approximation 
will  be  acceptable  and  linear  Kalman  Filter  theory  can  be  q)plied.  In  the  Extended 
Kalman  Filter,  the  plant  and  full  nonlinear  equations  of  state  are  used,  but  the  Kalman 
filter  gain  is  calculated  using  a  linear  approximation  of  the  model.  The  general  equations 
of  state  are  given  by: 

x(k  +  1)  =  f(x(k),  w(k),  k)  (28) 

Where:  x(k)  =  state  of  the  system,. 

w(k)=  the  plant  driving  noise. 

The  function  f(«)  can  be  a  nonlinear  function  of  the  states,  the  noise,  or  k. 
The  general  equations  for  the  observations  are  given  by: 

z(k)  =  h(x(k),v(k),k)  (29) 

Where:  z(k)  =  the  measurements  of  the  system,  and 

v(k)  =  the  measurement  noise  of  the  system. 
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The  function  h(*)  can  be  a  nonlinear  function  of  the  states,  the  measurement  noise,  or  k. 
The  measurement  noise  and  the  plant  driving  noise  are  assumed  to  be  zero  mean  white 

noise. 

The  extended  Kalman  filter  is  given  by  equations  (30)  through  (32).  These  equations 
utilize  the  nonlinear  relationships  to  predict  the  states  and  observations  of  the  the  system. 
The  smoothed  estimates  of  the  states  are  calculated  using  the  Kalman  gain  calculated 
with  a  linearized  version  of  the  Kalman  gain  equation. 

x(k+llk)  =  f(x(klk),k)  (30) 

z(k+llk)  =  h(x(k+llk),k)  (31) 

x(k+llk+l)  =  x(k+llk)+G(k+l)(z(k+l)-z(k+llk))  (32) 

Where  x(klk)  =  the  current  estimate  of  the  states, 

x(k+llk)  =  the  predicted  estimate  of  the  states, 
z(k  +  1)  =  the  observed  measurements, 

G(k  +1)  =  the  Kalman  gain, 

w(k) ,  the  plant  driving  noise,  is  assumed  to  be  zero  mean. 
v(k) ,  the  measurement  noise,  is  assumed  to  be  zero  mean. 

Equation  (30)  is  used  to  predict  the  next  state  of  the  system,  given  the  current  state  of 
the  system  and  the  current  input.  Since  the  plant  driving  noise  is  assumed  to  have  zero 
mean,  it  is  assumed  that  the  plant  driving  noise  does  not  contribute  the  expected  value  of 
the  next  state.  The  next  observation  is  predicted  using  equation  (3 1 ).  This  equation  uses 
the  predicted  state  of  the  system.  The  measurement  noise  is  assumed  to  have  zero  mean, 
and  does  not  contribute  to  the  the  expected  value  of  the  observation.  The  smoothed 
estimates  of  the  states  of  the  system  are  found  with  equation  (32).  This  equation  uses  the 
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Kalman  gain  ^^ch  is  calculated  using  equations  (33)  through  (35).  These  equations  use 
first  order  linear  approximations  of  state  prediction  and  observation  equations.  The 
covariance  prediction  equations  and  the  Kalman  gain  equations  are  given  as: 


P(k+llk)  =  6P(k!k)6'f  +  Q 

(33) 

P(k  + 1  Ik  + 1)  =  [l  -  G(k  +  l)H]p(k  + 1  Ik) 

(34) 

G(k  +  1)  =  P(k  +  1  lk)HT[HP(k  +  1  lk)HT  +  r]'' 

(35) 

Where  P(klk)  =  the  current  estimate  of  the  covariance  of  the  estimation  error, 

P(k+llk)  =  the  predicted  covariance,  given  the  current  value, 

P(k  +  1  Ik  +  1 )  =  the  predicted  covariance  given  the  current  observation, 
Q  =  the  covariance  of  the  plant  driving  noise, 

R  =  the  covariance  of  the  measurement  noise, 

4)  =  ^(x(k),u(k),w(k),k) 

x(k>sx(1clk),  u(k>=ii(k),  w(k>=0 

^  ^  gh(x(k),v(k),k) 

x(k>=x(klk).  v(k)=0 


To  initialize  the  Kalman  filter,  initial  estimates  for  the  states  and  the  covariance 
matrix  must  be  provided.  The  values  chosen  are  generally: 

x(OlO)  =  E[x(0)],  and  (36) 

P(0l0)  =  E[{x(0)-x(0l0)}{x(0)-x(0l0)}T]  =  cov[x(0)]  (37) 

The  Kalman  filter  provides  optimal  performance,  and  guarantees  convergence  and 
stability,  but  the  extended  Kalman  filter  does  not.  In  many  applications  the  extended 
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Kalman  filter  provides  accurate  estimates  of  the  states,  but  may  not  guarantee 
convergence  and  stability.  Generally  the  convergence  and  stability  are  highly  dependent 
on  the  values  chosen  for  the  initial  estimates  of  the  states,  and  covariance.  In  any 
application  of  the  extended  Kalman  filter,  the  convergence  and  stability  of  the  solutions 
must  be  thoroughly  explored  for  as  wide  a  range  of  initial  conditions  as  is  possible. 

B.  ERROR  ELLIPSOIDS 

The  error  ellipsoid  is  a  means  of  visualizing  the  error  in  the  estimate  of  the  states  of 
the  Kalman  filter.  The  uncertainty  in  the  estimate  is  expressed  in  the  error  covarianc';  P. 
The  error  covariance  matrix  P  is  the  expected  value  of  the  covariance  of  the  error  in  the 
states. 

P  =  E[(x-xKx-x)'^]  (38) 

Where  x  =  the  state  estimate  of  the  Kalman  filter, 

X  -  the  mean  of  the  states  of  the  Kalman  filter 

If  the  noise  entering  the  Kalman  filter  is  Gaussian  then  the  estimates  of  the  states  of 
the  Kalman  filter  also  have  a  Gaussian  distribution.  This  is  because  a  linear 
transformation  of  a  set  of  Gaussian  random  variables  will  result  in  another  set  of 
Gaussian  random  variables. 
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The  equation  for  the  probability  density  function  for  a  joint  Gaussian  random 


variable  is; 


'  (2n)’«id,.;r  ““ 


Where  x  =  the  vector  of  N  random  variables, 

X  =  the  expected  value  of  x,  E[x],  and 
P  =  the  error  covariance  matrix  of  X. 

The  error  covariance  matrix  is  given  by  the  following: 


Pll 

P«2 

•  •  •  Pin 

P21 

P22 

•  .  •  P2N 

Pni 

PN2 

.  .  .  Pnn 

(40) 


Where  py  =  E[(xi  -  Xi)(Xj  -  xj)] 

A  surface  in  N  space  can  be  found  \s1iere  the  probability  of  x  will  be  a  constant.  This 
surface  exists  over  the  values  of  x  for  which  the  argument  (x  -  x)^P“‘  (x  -  x)  is  a  constant. 
The  total  probability  that  the  values  of  x  are  lie  on  the  surface  defined  is  found  by 


integrating  the  probability  density  function  over  the  entire  surface.  These  surfaces  of 


constant  probability  are  called  error  ellipsoids  of  constant  probability.  If  the  argument 
(x  -  x)^P~‘ (x  -  x)  is  set  equal  to  the  constant  c^  then  the  probabilities  that  the  value  of  x 
will  lie  within  the  ellipsoid  formed  in  N  space  is  given  in  Table  1 . 
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c 

1 

2 

3 

N 

1 

0.683 

0.955 

0.997 

2 

0.394 

0.865 

0.989 

3 

0.200 

0.739 

0.971 

Table  1.  PROBABILITIES  FOR  ERROR  ELLIPSOIDS 


In  the  Kalman  filter  application  pursued  here  the  estimates  of  the  location  of  the 
emitter  are  joint  Gaussian  random  variables.  The  mean  of  the  steady  state  estimate  of  the 
location  of  the  emitter  is  considered  the  true  location  of  the  emitter,  and  the  difference 
between  the  estimate  of  the  location  of  the  emitter  and  the  mean  of  the  steady  state 
estimates  is  the  error.  This  error  term  is  used  to  form  the  error  covariance  in 
equation  (38),  Redefining  the  error  in  the  estimate  of  the  location  of  the  emitter,  the 
argument  of  the  joint  probability  density  function  is: 

c2=(xTp-'x)  (41) 


,  the  error  in  the  estimate  of  the  emitter  location, 

P  =  the  error  covariance  matrix 
c  =  a  constant. 

When  the  argument  in  equation  (41)  is  expanded  and  the  inverse  of  the  P  matrix  is 
found,  the  following  equation  results: 


Pll  Pl2 
Pl2  P22 


,  and 


Where 


X  = 


xt 

?h-xt 

. 

^-yt 

c^  = 


1 


P11P22  -P12 


-[ «  yt  ] 


'  P22  -P12 1 

1 - 

IX 

1 — 

1 

y 

— 1 

— j 

(42) 
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This  further  reduces  to; 


P22^^  -2Pu?hyt  +Pn^^ 

PnP22-PL 


= 


(43) 


To  facilitate  plotting  this  equation,  the  y  estimate  of  the  location  of  the  emitter  is  solved 

for  in  terms  of  the  x  coordinates  of  the  emitter.  Equation  (43)  is  rewritten  in  the  form  of  a 
quadratic  equation  Ay^  +  By  +  C  =  0  wdiere 


A  = 


11 


PllP22-Pt2 


B  = 


-2P,2Xt 


C - 


P11P22-P12  P11P22-P12 


The  quadratic  formula  is  used  to  calculate  the  value  of  ytin  terms  of  xt  and  the  elements 
of  the  error  covariance  matrix.  The  following  equation  results: 


yt  = 


Pl2Xt^  |(P||P22-Pj2)^  _2 


11 


-(PllC^-XtO 


(44) 


11 


The  maximum  and  minimum  values  of  the  error  xtare  found  from  the  locations  where 
xt^  =  Pnc^ .  If  xt^  >  Pnc^  a  real  solution  to  equation  (44)  will  not  exist. 
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V.  SIMULATION  RESULTS 
A.  THE  BURST  TDOA  KALMAN  FILTER 


1.  Extended  Kalman  Filter  Equations 

An  extended  Kalman  filter  algorithm  is  developed  to  estimate  of  the  emitter 
location  from  burst  TDOA  observations.  Simulated  data  is  used  to  test  the  accuracy  and 
stability  of  the  algorithm.  The  TDOA  measurement  noise  is  assumed  to  be  white 

gaussian.  The  following  assumptions  are  made  to  simplify  the  problem: 

1 .  The  emitter  scans  in  azimuth  at  a  constant  rate. 

2.  An  estimate  of  the  scan  rate  of  the  emitter  is  known  a  priori. 

3.  Some  method  is  used  to  calculate  the  TDOA  of  the  radar  signal  between  the  two 
receivers,  and  the  error  on  this  TDOA  observation  is  modeled  as  additive  white 
Gaussian  noise  with  variance  R. 

4.  The  emitter  is  stationary. 

5.  The  X  and  y  coordinates  of  the  emitter  and  receivers  are  measured  with  respect  to 
a  fixed  coordinate  reference. 

6.  The  exact  locations  of  the  receivers  are  known  for  all  points  in  time. 

These  assumptions  simplify  the  problem,  but  do  not  reduce  its  usefulness,  because  the 

algorithm  can  easily  be  modified  to  account  for  these  assumptions. 

The  estimated  state  and  observations  are  calculated  from  the  following: 

x(k+l)  =  x(k)=  (45) 


z.b(k)  = 


J _ (xt  -  xa)(yt  -  yb)  -  (yt  -  ya)(xt  -  xb) 

[[(xt  -  xa)^  +  (yt  -  ya)^][(xt  -  xb)^  +  (yt  -  yb)^]] 


+  V|*(k) 


(46) 
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+  VkOc)  (47) 


^  J _ (xt  - xaXyt - yc) - (yt -  ya)(xt  -  xc) 

[[(xt  -  xa)^  +  (yt  -  ya)^][(xt  -  xc)^  +  (yt  -  yc)^]] 

i  1  (xt-xb)(yt-yc)-(yt-ybXxt-xc)  , 

Zbc(K)  =  —  - ; - ;; - ; - — —  +  Vbc(k)  (48) 

[[(xt  -  xhy  +  (yt  -  yb)^][(xt  -  xc)^  +  (yt  -  yc)^]]  J 

Where  xt,yt  =  the  estimate  of  the  emitter's  x  and  y  coordinates, 

SR  =  the  scan  rate  of  the  emitter, 

v„  =  the  measurement  noise  present  in  the  TDOA  observation 
between  two  receivers, 

xa,  ya  =  the  X  and  y  coordinates  of  receiver  A, 

xb,  yb  =  the  X  and  y  coordinates  of  receiver  B, 

xc,  yc  =  the  X  and  y  coordinates  of  receiver  C, 

xt,  yt,  xa,  ya,  xb,  yb,  xc,  and  yc  are  all  functions  of  k. 

Since  the  emitter  is  stationary,  the  state  transition  matrix  in  equation  (45)  is  the 
identity  matrix.  Therefore,  the  predicted  estimate  of  the  emitter's  location  is  the  same  as 
the  current  estimate  of  the  location.  Since  the  coordinates  of  receivers  A,  B,  and  C  are 
known  at  all  points  in  time  and  the  scan  rate  is  known  a  priori,  the  predicted  TDOA 
observations  are  only  functions  of  the  predicted  states. 

The  Kalman  filter  is  implemented  as  three  separate  Kalman  filters  linked 
together.  Each  receiver  processes  the  received  observations  separately.  They  calculate 
Kalman  gains  and  an  error  covariance  matrix  from  its  observations,  and  then  use  these  to 
calculate  an  estimate  of  the  location  of  the  emitter.  This  estimate  is  then  passed  to  the 
next  receiver.  This  algorithm  is  used  because  it  facilitates  implementation  in  hardware. 
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The  equations  necessary  for  the  Extended  Kalman  Filter  implementation  are: 


P(k+llk)  =  P(klk)  +  Q  (49) 

P(k  + 1  Ik  + 1)  =  [l  -  G«(k  +  l)H«]p(k  + 1  Ik)  (50) 

Gxx(k  + 1)  =  P(k  + 1  lk)HT  [HxxP(k  + 1  \k)Hl,  +  V**  ]"'  (51) 

Where  P(k+l|k)  =  the  predicted  error  covariance, 

P(k|k)  =  the  current  error  covariance, 

Q  =  the  covariance  of  the  state  disturbances, 

G„(k+1)  =  the  Kalman  Gain  calculated  for  the  TDO A  observation 
between  tv;o  of  the  receivers, 

V„  =  the  covariance  of  the  measurement  noise  for  the  TDOA 
observation  between  two  of  the  receivers, 

„  ^  r  ai„oc)  3ho(k)  I 

"w  1_  3yt  J 

The  algebraic  equations  that  make  up  the  partial  derivatives  in  H„  are  extensive  and  are 
omitted  here  fo»-  cb  The  derivation  of  is  included  in  Appendix  D. 

The  flow  chart  for  the  Kalman  filtering  algorithm  is  shown  in  Figure  13.  This 
flowchart  provides  the  basic  logic  flow  of  the  Kalman  filtering  algorithm.  In  an  actual 
implementation,  TDOA  observations  are  available,  and  only  the  portion  contained  within 
the  dotted  outline  is  required.  The  additional  portion  of  the  algorithm  generates  the 
observations  used  to  test  the  algorithm. 
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2.  Burst  TDOA  Simulations 


The  following  scenario  is  used  to  test  the  burst  TDOA  Kalman  filter  algorithm. 
Three  receivers  are  used  and  all  three  of  the  available  TDOA  observations  are  used  to 
estimate  the  location  of  the  emitter.  The  initial  locations  of  the  receivers  and  their 
velocities  and  bearings  are  shown  in  Figure  14. 


> 

J 

< 

Receiver  A 
(0.0)  km 

^  V  =  30  m/sec 

Bearing  90  deg. 

w  V 

Receiver  C 
(-20.-20)  km 

V  =  42  m/sec 

Bearing  -135  deg. 

/ 

\ 

Receiver  B 
(20.-2)  Km 

V  *  42  m/sec 

Bearing  -45  deg 

Figure  14.  Initial  location  and  velocity  of  receivers 


The  emitter  locations  are  shown  in  Figure  15. 
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Although  the  aigorithm  has  the  capability  to  account  for  different  signal  strengths 
and  resulting  error  variances,  all  of  the  receivers  are  considered  identical  with  identical 
signal  to  noise  ratios.  An  8  MHz  sampling  rate  is  assumed.  The  algorithms  developed  in 
Appendix  A  are  used  to  calculate  the  variance  of  the  error  present  in  the  TDOA 
observations  for  the  signal  to  noise  ratio  and  the  'ompling  rates  chosen.  In  the  tests  of 
Kalman  filtering  algorithm,  this  error  is  modeled  as  zero  mean  white  Gaussian  noise,  that 
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is  added  to  the  true  TDOA  observations.  The  variance  of  the  TDOA  observation  noise  is 


5.660  X  10  ’  for  all  of  the  scenarios. 

For  all  of  the  scenarios  the  a  priori  estimate  of  the  emitter  location  is  chosen  as 
(0,10)  kilometers.  The  filter  parameters  Q,  the  covariance  of  the  state  excitation,  and  Po, 
the  initial  error  covariance,  are  varied  until  the  filter  converges  well.  All  of  the  scenarios 
are  run  for  twenty  five  TDOA  observations.  Given  the  scan  rate  of  the  emitter,  the 
real-time  length  of  the  simulation  is  twenty  five  seconds. 

Three  plots  are  presented  for  each  scenario.  The  first  plot  is  an  X-Y  plot  of  the 
estimates  of  the  coordinates  of  the  emitter.  This  plot  demonstrates  the  track  of  the  final 
estimate  of  the  emitter  location  after  all  three  TDOA  observations  are  processed  for  each 
time  step.  In  this  plot  the  loci  of  constant  TDOA  that  correspond  to  the  final  steady  state 
TDOA  observations  are  also  plotted.  Some  of  the  loci  do  not  intersect  exactly  at  the 
emitter.  The  equations  used  to  plot  the  loci  assume  that  the  distance  to  the  emitter  is 
much  greater  than  the  separation  of  the  receivers,  and  in  some  of  the  scenarios  this 
assumption  is  not  valid.  The  loci  visually  illustrate  the  relationship  between  the  TDOA 
observations  tmd  give  a  visual  indication  of  the  orthogonality  of  the  observations.  The 
3o  error  ellipsoids  are  also  plotted  for  the  1st,  8th,  15th,  and  22nd  estimates  of  the  emitter 
location.  These  ellipsoids  provide  a  visual  indication  of  the  accuracy  of  the  location 
estimate. 


44 


The  second  plot  is  the  estimate  of  the  x  and  y  coordinates  of  the  emitter  vs.  the 
number  of  TDOA  observations.  This  plot  gives  a  visual  indication  of  how  rapidly  the 
estimates  of  the  emitter  location  converge  on  the  true  coordinates  of  the  emitter. 

The  third  plot  is  an  X-Y  plot  that  shows  the  true  emitter  coordinates,  the  loci  of 
constant  TDOA,  and  the  estimates  of  the  emitter  coordinates  in  the  vicinity  of  the  mean 
steady  state  estimate  of  the  emitter  coordinates.  The  3o  error  ellipsoid  is  plotted  for  the 
steady  state  estimate,  and  the  length  of  the  major  and  minor  axis  of  the  error  ellipsoid 
give  an  indication  of  the  maximum  error  in  the  estimate. 
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a.  Scenario  U1 


In  the  first  scenario  the  performance  of  the  burst  TDOA  filter  is  tested  for  an 
emitter  at  a  range  of  SOO  kilometers  from  the  origin  and  a  bearing  of  90  degrees  from  the 
X  coordinate  axis.  The  initial  error  covariance  matrix,  Po,  and  the  variance  of  the  state 
excitation,  Q,  are  varied  until  the  filter  converges  to  a  steady  state  value  within  ten 
observations.  The  plots  presented  in  figures  16  through  18  are  representative  of  the 
results  obtained.  The  filter  converges  very  well,  and  as  shown  in  figures  16  and  18  the 
estimate  of  the  location  of  the  emitter  converges  to  the  true  location  of  the  emitter.  The 


values  that  gave  the  best  convergence  were: 


Po  = 


1.0  0 
0  1.0 


4.0  0 
0  4.0 


X  10"m^ 


xl0“m^and  Q  = 

As  shown  in  Figure  16,  the  3a  error  ellipsoids  grow  larger  as  the  estimate  of 
the  location  of  the  emitter  grows  closer  to  the  actual  location.  Because  the  distance  to  the 
emitter  is  large  in  relation  to  the  separation  of  the  receivers,  the  TDOA  observations  are 
very  close  to  each  other.  The  estimated  location  of  the  emitter  does  not  approach  the  true 
location  until  the  error  covariance  matrix  and  subsequently  the  Kalman  gains  grow  large 
enough  to  amplify  the  very  small  differences  in  the  TDOA  and  drive  the  states  of  the 
filter.  A  high  value  for  Q  is  required  to  drive  the  error  covariance  matrix  high  enough.  If 
a  smaller  value  of  Q  is  used,  the  filter  reached  a  steady  state  value  far  from  the  true 
location  of  the  emitter. 
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Figure  16.  Scenario  #1  Estimates  of  the  location  of  the  emitter 
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Figure  18.  Scenario  #1  Close-up  of  steady  state  estimate  of  emitter  location 
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b.  Scenario  if 2 


In  the  second  scenario,  the  performance  of  the  burst  TDOA  filter  is  tested  at 
an  intermediate  range.  The  emitter  location  is  chosen  at  a  range  of  150  kilometers  from 
the  origin  and  at  a  bearing  of  30  degrees  from  the  x  coordinate  axis.  The  coordinates  of 
the  emitter  are  (130, 75)  kilometers.  The  initial  error  covariance  matrix,  Po,  and  the 
variance  of  the  state  excitation  term,  Q,  are  varied  until  the  filter  converges  to  a  the 
correct  location  within  approximately  ten  observations.  The  plots  presented  in  figures  19 
through  21  are  representative  of  the  results  obtained.  The  values  that  give  the  best  results 
are: 


Po 


1.0  0 
0  1.0 


5.0  0 
0  5.0 


X  lO^m^ 


xl0^°m^and  Q  = 

The  filter  performed  very  well.  The  orthogonality  of  the  TDOA  observations 
in  the  vicinity  of  the  emitter  is  high,  and  the  estimate  of  the  emitter  location  converged  to 
the  proper  location  quickly  and  smoothly.  The  final  steady  state  estimate  is  stable  and  no 
numerical  problems  are  present.  The  3a  error  ellipsoids  are  not  visible  in  figure  19 
because  of  the  scale  of  the  figure. 
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y  coordinate  (kilometers) 


Figure  19.  Scenario  #2  Estimates  of  emitter  location 
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Figure  20.  Scenario  #2  Estimates  of  emitter  coordinates  vs  observations 
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Figure  21.  Close-up  of  the  steady  state  estimate  of  emitter  location 


c.  Scenario  #5 


In  the  third  scenario  the  performance  of  the  burst  TDOA  filter  is  tested  for  an 
emitter  at  close  range.  The  emitter  location  is  chosen  at  a  range  of  30  kilom^ers  fi-om  the 
origin  and  at  a  bearing  of  60  degrees  fi-om  the  x  coordinate  axis.  The  initial  error 
covariance  matrix,  Po,  and  the  variance  of  the  state  excitation  term,  Q,  are  varied  imtil  the 
filter  converged  to  the  correct  location  of  the  emitter  within  approximately  ten 
observations.  The  plots  presented  in  figures  22  through  24  are  representative  of  the 


results  obtained.  The  values  that  give  the  best  results  are; 


1.0  0 
0  1.0 


X  10‘*  m^  and  Q  = 


40  0 
0  40 


Although  the  loci  of  constant  TDOA  plotted  in  figure  22  do  not  intersect  at 


the  emitter,  they  indicate  that  the  orthogonality  of  the  TDOA  observations  is  high  and 


accurate  results  are  expected.  Figure  23  shows  that  the  estimates  of  the  emitter  location 


converged  to  the  correct  location  of  the  emitter  rapidly  and  smoothly.  Because  of  the 


large  scale  of  the  plot,  the  3a  error  ellipsoids  are  not  visible  in  figure  22. 
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X  coordinate  (kilometers) 


Figure  22.  Scenario  #3  Estimates  of  emitter  location 


Figure  23.  Scenario  #3  Estimates  of  emitter  coordinates  vs  observations 


Figure  24.  Close-up  of  the  steady  state  estimate  of  emitter  location 
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3.  The  Burst  TDOA  Filter  Results 


The  simulations  presented  here  in  no  way  test  all  the  possible  implementations  of 
this  filtering  algorithm.  Further  evaluation  of  the  algorithm  is  necessary  to  fully  evaluate 
its  performance  and  capabilities  given  the  wide  variety  of  possible  emitter  locations, 
receiver  configurations,  and  emitter  characteristics.  From  the  an<ilysis  conducted  some 
general  observations  can  be  drawn. 

For  the  burst  TDOA  Kalman  filter  the  orthogonality  among  the  TDOA 
observations  is  generally  very  good.  For  the  three  randomly  chosen  receiver  locations, 
the  orthogonality  of  the  TDOA  observations  are  not  heavily  dependent  on  the  location  of 
the  emitter.  As  the  locations  of  the  receivers  remained  constant  and  the  bearing  to  the 
emitter  is  varied  throughout  a  90  degree  arc,  and  there  does  not  appear  to  be  any 
particular  bearings  where  the  orthogonality  of  all  three  observations  would  decrease  to 
the  point  where  the  performance  of  the  filter  degrades.  For  all  bearings  and  ranges 
examined  at  least  two  of  the  TDOA  observations  had  good  orthogonality. 

The  accuracy  of  the  filter  is  heavily  dependent  on  the  values  chosen  for  Q,  the 
covariance  of  the  state  excitation.  If  this  value  is  too  low,  the  filter  reaches  a  steady  state 
value,  but  the  steady  state  value  is  not  the  correct  location  of  the  emitter.  As  Q  is 
increased,  the  steady  state  estimate  of  the  location  of  the  emitter  reaches  a  stable  point 
around  the  true  emitter  location.  The  magnitude  of  Q  directly  affects  the  magnitude  of 
the  error  covariance  P,  and  as  Q  is  increased  the  error  covariance  and  the  3a  error 
ellipsoids  for  the  steady  state  estimate  increase  in  size. 
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The  filter  performs  well  for  the  receiver  configurations  and  emitter  locations 
tested.  Even  at  the  maximum  range,  where  the  receivers  are  spaced  closely  in  comparison 
to  the  distance  from  the  emitter,  the  error  present  in  the  final  estimate  of  the  location  of 
the  emitter  is  reasonable.  For  scenario  #1  the  length  of  the  3a  error  ellipsoid,  which 
represents  effectively  the  maximum  error  in  the  estimate  given  the  noise  present  in  the 
observations,  is  roughly  30  kilometers  in  bearing  direction  and  10  kilometers  in  range. 
These  represent  errors  of  about  3  degrees  in  bearing  and  2  percent  in  range.  In  scenarios 
#2  and  #3,  the  distance  to  the  emitter  decreases  relative  to  the  separation  of  the  receivers 
and  the  maximum  error  decreases  to  about  0.5  and  0.20  degrees  in  bearing,  and  2  percent, 
and  0.3  percent  in  range  respectively. 

These  results  are  accomplished  filtering  the  TDOA  observations  one  at  a  time. 
This  filtering  technique  tends  to  skew  the  error  ellipsoids  to  align  them  vrith  the  loci  of 
constant  TDOA  that  corresponded  to  the  last  TDOA  observation.  This  places  the  major 
axis  of  the  error  ellipsoid  along  this  loci  and  tends  to  increase  the  error  along  that  loci.  If 
all  three  of  the  TDOA  observations  are  processed  simultaneously  as  the  filter  reaches  a 
steady  state  solution,  the  error  covariance  and  thus  the  error  ellipsoids  and  maximum 
error  present  in  the  estimate  can  be  decreased  further. 
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B.  SINGLE  PULSE  TDOA  KALMAN  FILTER 


1.  Extended  Kalman  Filter  Equations 

An  extended  Kalman  filter  algorithm  is  developed  to  estimate  of  the  emitter 
location  fix>m  pulse  TDOA  observations.  Simulated  data  is  used  to  test  the  accuracy  and 

stability  of  the  algorithm.  The  following  assumptions  are  made  to  simplify  the  problem: 

1 .  Two  receiver  platforms  are  used,  and  each  receiver  platform  is  equipped  with  two 
sensors  to  detect  the  pulse  TOA. 

2.  The  error  present  in  the  pulse  TDOA  observation  is  modeled  as  zero  mean 
additive  white  Gaussian  noise  with  variance  R. 

3.  The  emitter  is  stationary. 

4.  The  X  and  y  coordinates  of  the  emitter  and  receivers  are  measured  with  respect  to 
a  fixed  coordinate  system,  and  the  coordinates  of  the  sensors  are  known  exactly. 

5.  The  PRI  is  long  enough  that  ambiguity  is  not  present  in  the  pulse  TDOA 
observations. 

These  assumptions  simplify  the  problem,  but  do  not  reduce  its  usefulness,  because  the 
algorithm  developed  can  easily  be  modified  to  account  for  these  assumptions. 

The  estimated  state  and  predicted  observation  are  calculated  fi'om  the  following: 

x(k+l)  =  x(k)=  (52) 


z(k+llk)  = 

The  Kalman  filter  equations  are: 

x(k+llk+l)  =  x(k+llk)  +  G(k+l)[z(k+l)-z(k+llk)]  (54) 

P(k+llk)  =  P(klk)  +  Q  (55) 


(xt(k)  -  xa(k))^  +  (yt(k)  -  ya(k))^  1  -  [  (xt(k)  -  xb(k))^  +  (yt(k)  -  yb(k))^ 
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(56) 


P(k  +  1  Ik  1)  =  [l  -  G(k  +  l)H(k)]p(k  + 1  Ik) 

G(k  +  1)  =  P(k  +  1  lk)H(k)T[H(k)P(k  +  1  lk)H(k)T  +  r]'‘  (57) 

Where  xt(k)  ^(k  =  the  estimate  of  the  x  and  y  coordinates  of  the  emitter, 

z(k+l)  *=  the  TDOA  observations  for  receiver  platforms  1  and  2, 
z(k  + 1  Ik)  =  the  predicted  TDOA  based  upon  the  estimates  of  emitter 
location. 

xa,  ya  =  the  X  and  y  coordinates  of  the  .sensor  A  of  the  receiver  platform, 

xb,  yb  =  the  X  and  y  coordinates  of  the  sensor  B  of  the  receiver  platform, 
G(k+1)  =  the  Kalman  gains, 

P(k|k)  =  the  error  covariance, 

H(k)  =  the  gradient  of  the  observation  equation, 

Q  =  the  covariance  of  the  plant  noise, 

R  =  the  covariance  of  the  measurement  noise, 

c  =  the  speed  of  light, 

xa,  xb,  ya,  yb  are  all  functions  of  k. 

The  gradient  of  the  observation  equation  H(k)  is  calculated  from  the  following  equations. 


TT/n.x  r  lf(xt(k>-xa)  (xt(k>-xb)'l  lf(yt(k)-ya)  (yt(k>-yb)'\ 

cV  Ra(k)  Rb(k)  J  Ra(k)  Rbflc)  J 

(58) 

Where, 

Ra(k)  =  [(5rt(k)  -  xa)^  +  (^(k)  --  ya)^  J  ^  ,  and 

(59) 

Rb(k)  =  [(xt(k)  -  xb)^  +  (>h(k)  -  yb)^  ] 

(60) 

The  flowchart  for  the  pulse  TDOA  Kalman  filter  algorithm  is  identical  to  the 
flowchart  used  for  the  burst  TDOA  Kalman  filter  shown  in  Figure  13.  The  simulation 
calculates  the  pulse  TDOA  observation  using  the  known  locations  of  the  emitter  and 
sensors.  Zero  mean  white  Gaussian  noise  is  ar*  Jed  to  the  calculated  observations  and  the 
filter  estimates  the  location  of  the  emitter  from  these  noisy  observations.  The  variance  of 
the  measurement  noise  is  calculated  from  the  emitter  parameters  and  an  assumed  signal  to 
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noise  ratio.  The  filtering  algorithm  is  tested  for  a  niunber  of  emitter  locations.  In  each 
scenario  the  initial  error  covariance  matrix  and  the  variance  of  the  plant  excitation  noise 
are  varied  until  the  filter  converges  to  the  emitter’s  true  coordinates 

2.  Pulse  TDOA  Kalman  Filter  Simulations 

Three  scenarios  are  used  to  test  the  performance  of  the  pulse  TDOA  Kalman 
filter.  The  receiver  locations  remain  the  same  for  all  of  the  scenarios  and  are  shown  in 
Figure  25.  The  sensors  are  stationary  for  all  of  the  scenarios.  The  sensors  are  separated 
by  500  meters,  and  the  receiver  platform  separation  is  20  kilometers. 
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Figure  25.  Receiver  platform  and  sensor  configuration 
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The  emitter  locations  chosen  to  test  the  pulse  TDOA  algorithm  are  shown  in  Figure  26. 


Figure  26.  Emitter  locations  for  the  single  pulse  TDOA  Kalman  filter 


The  emitter  has  the  following  characteristics.  These  characteristics  are  typical  for  long 


range  search  radar. 

Beam  width 
Scan  Rate 
PRF 

Pulsewidth 

Peak  Signal  to  Noise  Ratio 


2.0  deg. 

360  deg/sec. 
2000  Hz 
1 .0  microsecond 
15dB 


All  of  the  sensors  are  considered  identical  with  identical  peak  signal  to  noise 
ratios.  An  8  MHz  sampling  rate  is  assumed.  The  algorithms  developed  in  Appendix  A 
are  used  to  calculate  the  variance  of  the  error  present  in  the  TDOA  observations  for  the 
sampling  rate  and  signal  to  noise  ratio  chosen.  In  the  tests  of  the  Kalman  filtering 
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algorithm,  this  error  is  modeled  as  zero  m  -isn  white  Gaussian  noise,  and  is  added  to  the 
true  TDOA  observations.  The  variance  of  the  TDOA  observation  noise  is  1.694  x  10  '* 
for  all  of  the  scenarios. 

For  all  of  the  scenarios  the  a  priori  estimate  of  the  emitter's  location  is 
(10,5)  kilometers.  This  a  priori  estimate  provided  good  convergence  for  a  wide  range  of 
emitter  locations.  The  variance  of  the  plant  excitation  noise  Q,  and  the  initial  error 
covariance  are  varied  until  the  filter  converges  to  the  actual  coordinates  of  the  emitter. 
All  three  of  the  scenarios  are  run  for  sixty  TDOA  observations.  For  the  emitter 
characteristics  chosen,  this  represents  a  real  run  time  of  30  milliseconds. 

Three  plots  are  presented  for  each  of  the  scenarios  examined.  The  first  plot  is  an 
X- Y  plot  of  the  estimates  of  the  coordinates  of  the  emitter.  This  plot  demonstrates  the 
track  of  the  final  estimate  of  the  emitter  location  after  both  TDOA  observatioas  are 
processed  at  each  time  step.  The  loci  of  constant  TDOA  that  correspond  to  the  final 
steady  state  TDOA  observations  are  also  plotted.  These  loci  illustrate  the  relationship 
between  the  TDOA  observauons  and  give  a  visual  indication  of  the  orthogonality  of  the 
observations.  Also,  the  3a  error  ellipsoids  are  plotted  for  the  1st,  20th,  39th,  and  58th 
estimates  of  the  emitter  location.  These  ellipsoids  provide  a  visual  indication  of  the 
accuracy  of  the  location  estimate. 

The  second  plot  is  the  estimate  of  the  x  and  y  coordinates  of  the  emitter  vs.  the 
number  of  TDOA  observations.  This  plot  gives  a  visual  indication  of  how  rapidly  the 
states  of  the  filter  converge  to  the  true  coordinates  of  the  emitter. 
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The  third  plot  is  an  X-Y  plot  of  true  emitter  coordinates,  the  loci  of  constant 
TDOA,  and  the  estimates  of  the  emitter  coordinates  in  the  vicinity  of  the  mean  steady 
state  estimate  of  the  emitter  coordinates.  The  3a  error  ellipsoid  is  plotted  for  the  steady 
state  estimate.  The  length  of  the  major  and  minor  axis  of  the  error  ellipsoid  give  an 
indication  of  the  maximum  error  in  the  estimate. 
a.  Scenario  U1 

In  the  first  scenario  the  performance  of  the  pulse  TDOA  filter  is  tested  for  an 
emitter  located  at  a  range  of  500  kilometers  fiom  the  origin  at  a  bearing  of  90  oegrees 
fi-om  the  X  coordinate  axis.  The  initial  error  covariance  matrix,  Po,  and  the  variance  of 
the  state  excitation,  Q,  are  varied  until  the  filter  converges  to  a  steady  state  value.  The 
plots  presented  in  Figures  27  through  29  are  representative  of  the  results.  In  each 
simulation,  a  total  of  sixty  observations  are  filtered.  The  values  that  give  the  best  results 
are: 


Po  = 


5.0  0 
0  5.0 


xlO^^m^and  Q  = 


1.0 

0 


0 

1.0 


X  10*  m^ 
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Figure  27.  Scenario#!  Estimates  of  emitter  location 


The  filter  performs  well,  but  the  error  in  the  steady  state  estimate  is  high,  and 
the  estimate  of  the  y  coordinate  of  the  emitter,  shown  in  Figure  28,  converges  very 
slowly.  Because  the  distance  to  the  emitter  is  large  in  relation  to  the  separation  of  the 
receivers,  the  TDOA  observations  are  very  close  to  each  other.  The  estimated  location  of 
the  emitter  does  not  approach  the  true  location  of  the  emitter  until  the  error  covariance 
matrix  and  subsequently  the  Kalman  gains  grow  large  enough  to  amplify  the  very  small 
differences  in  the  TDOA  to  drive  the  states  of  the  filter  near  the  true  location.  A  large 
value  for  Q  is  required  to  drive  the  error  covariance  matrix  high  enough.  If  a  smaller 
value  of  Q  is  used,  the  filter  reaches  a  steady  state  far  from  the  true  location  of  the 
emitter.  Because  of  the  magnitude  of  the  error  covariance  matrix,  the  3o  error  ellipsoids 
grow  large.  The  major  axis  of  the  steady  state  error  ellipsoid  is  nearly  400  kilometers. 
Considering  that  the  distance  to  the  emitter  is  only  500  kilometers,  a  location  estimate 
with  that  much  possible  error  has  limited  usefulness. 

In  a  second  attempt,  the  separation  of  the  receiver  platforms  is  increased  to 
300  kilometers  to  determine  if  greater  separation  of  the  receiver  platforms  can  increase 
the  orthogonality  of  the  TDOA  observations  and  improve  the  performance  of  the  filter. 
The  separation  of  the  sensors  is  maintained  at  500  meters.  The  results  of  the  simulation 
are  shown  in  Figures  30  through  32.  The  value  chosen  for  Q  in  this  simulation  is 
x  10*’.  The  value  chosen  for  Po  is  the  same  as  the  previous  simulation. 
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Figure  31.  Scenario  #1  Estimates  of  emitter  coordinates  vs.  observations  for 

300  km  receiver  separation 


Figure  32.  Scenario  #1  Close-up  of  steady  state  estimate  of  emitter  location  for  300  km 

receiver  separation 
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As  shown  in  Figure  3 1  the  estimates  of  the  emitter  coordinates  approached  the 
true  coordinates  of  the  emitter  more  quickly  and  directly  than  the  previous  simulation. 

This  is  largely  due  to  the  greater  oithogonahty  of  the  TDOA  observations.  The  steady 
state  3a  error  eUipsoid  plotted  in  Figure  32  is  considerably  smaller  than  in  the  previous 
simulation  and  the  major  axis  of  the  error  eUipsoid  is  only  S 1  kilometers. 
b.  Scenario  U2 

In  the  second  scenario  the  performance  of  the  burst  TDOA  filter  is  examined 
with  the  emitter  located  at  an  intermediate  range.  An  emitter  is  located  at  a  range  of  1  SO 
kilometers  from  the  origin  and  a  bearing  of  30  degrees  fi^om  the  x  coordinate  axis.  The  x 
and  y  coordinates  of  the  emitter  are  (130,  75)  kilometers.  As  in  the  first  simulation,  the 
initial  error  covariance  matrix,  Po,  and  the  variance  of  the  state  excitation  term,  Q,  are 
varied  until  the  filter  converges  to  the  true  coordinates  of  the  emitter.  The  values  that  give 
the  best  convergence  are: 


Po  = 


1.0  0 
0  1.0 


xlO  m  ,  and  Q  = 


1.0  0 
0  1.0 


xl0’m^ 


The  plots  presented  in  Figures  33  through  35  are  representative  of  the  results 


obtained. 
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y  coordinate  (kilometers) 


Figure  33.  Scenario  #2  Estimates  of  emitter  location 
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Figure  34.  Scenario  #2  Estimates  of  emitter  coordinates  vs.  observations 


Figure  35.  Close-up  of  the  steady  state  estimate  of  emitter  location. 
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The  filter  performs  well,  but  the  error  in  the  steady  state  estimate  is  high,  and 
the  estimates  of  the  x  and  y  coordinates  of  the  emitter,  shown  in  Figure  34,  converge  very 
slowly.  The  distance  fi’om  the  receiver  platforms  to  the  emitter  is  much  less  than  in 
scenario  #1,  but  the  orthogonality  of  the  TDOA  observations  is  low  because  of  the 
location  of  the  emitter  in  relation  to  the  receivers.  As  shown  in  figures  33  and  35,  the  loci 
of  constant  TDOA  are  nearly  parallel  in  the  vicinity  of  the  emitter.  Although  the  filter 
converges  to  a  steady  state  value  close  to  the  actual  emitter  location,  the  error  covariance 
is  very  high.  The  major  axis  of  the  3o  error  ellipsoid  for  the  steady  state  estimate  is 
approximately  140  kilometers.  Considering  that  the  distance  to  the  emitter  is  only  150 
kilometers,  a  location  estimate  with  an  error  ellipsoid  this  large  has  limited  usefulness. 

As  shown  in  the  first  scenario,  greater  separation  of  the  receiver  platforms  can 
increase  the  orthogonality  of  the  TDOA  observations  and  improve  the  performance  of  the 
filter.  When  the  separation  of  the  receiver  platforms  is  increased  to  1 50  kilometers  and 
the  separation  of  the  sensors  maintained  at  500  meters,  the  performance  improves 
dramatically.  The  results  of  the  simulation  are  shown  in  figures  36  through  38.  The  final 

r  5.0  0  1  , 

value  chosen  for  Q  is  q  5  q  ^  value  chosen  for  Po  is  the  same  as  the 

previous  simulation. 
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y  coordinate  (kilometers) 


X  coordinate  (kilometers) 


Figure  36.  Scenario  #2  Estimates  of  the  location  of  the  emitter  with 
1 50  km  receiver  separation 
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Figure  38.  Close-up  of  the  steady  state  estimate  of  emitter  location  with  1 50  km  receiver 

separation 
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As  shown  in  Figure  37,  the  estimates  of  the  emitter  coordinates  approach  the 
true  values  more  quickly  and  directly  than  the  previous  simulation.  This  is  largely  due  to 


the  greater  orthogonality  of  the  TDOA  observations.  The  3o  error  ellipsoid  plotted  in 
Figure  38  is  considerably  smaller  than  in  the  previous  simulation  The  major  axis  of  the 
error  ellipsoid  is  only  16.7  kilometers, 
c  Scenario  #3 

In  the  third  scenario  the  performance  of  the  pulse  TDOA  filter  is  tested  for  an 
emitter  located  at  close  range  in  comparison  to  the  separation  of  the  receivers.  The 
emitter  location  is  located  at  a  range  of  30  kilometers  fi'om  the  origin  and  at  a  bearing  of 
60  degi  i-es  from  the  x  coordinate  axis.  The  x  and  y  coordinates  of  the  emitter  are  ( 1 5,  26) 
kilometers  As  in  the  previous  simulations,  the  initial  error  covariance  matrix  and  the 
vanance  of  the  state  excitation  term,  are  varied  until  the  filter  converges  to  the  true 
coordinates  of  the  emitter.  The  values  that  give  the  best  convergence  are: 
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The  plots  presented  in  Figures  39  through  41  are  representative  of  the  results 

obtained. 

As  shown  in  Figure  40  the  estimate  of  the  emitter  location  converges  to  the 
true  coordinates  of  the  emitter  quickly.  Increasing  the  value  of  Q  improves  the 
convergence  of  the  filter,  but  also  increases  the  error  present  in  the  final  estimate  and 
increases  the  size  of  the  error  ellipsoid  For  the  values  of  Q  chosen  the  major  axis  of  the 
3o  error  ellipsoid  is  only  1.4  kilometers 
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Figure  39.  Scenario  #3  Estimates  of  the  emitter  location. 
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Figure  40.  Scenario  #3  Estimates  of  emitter  coordinates  vs.  observations 


Figure  41.  Close-up  of  the  steady  state  estimate  of  the  emitter  location 
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3.  Pulse  TDOA  Filter  Results 


The  simulations  presented  here  do  not  test  all  the  possible  implementations  of 
this  filtering  algorithm.  Further  evaluation  of  the  algorithm  is  necessary  to  fully  evaluate 
its  performance  and  capabilities  given  the  infinite  variety  of  possible  emitter  locations, 
receiver  configurations,  and  emitter  characteristics.  From  the  analysis  conducted  some 
general  observations  can  be  drawn. 

The  pulse  TDOA  filter  accuracy  and  effectiveness  is  heavily  dependent  on  the 
orthogonality  of  the  TDOA  observations.  As  the  orthogonality  of  the  TDOA 
observations  decreased,  the  error  present  in  the  estimate  increased.  When  the  receiver 
platforms  are  closely  spaced  relative  to  the  distance  to  the  emitter,  the  loci  are  nearly 
parallel  to  each  other  and  the  TDOA  observations  have  low  orthogonality.  To  obtain 
good  filtering  results  the  spacing  between  the  receiver  platforms  has  to  be  increased. 

The  receiver  locations  chosen  for  the  single  pulse  TDOA  filtering  problem 
yielded  areas  where  the  orthogonality  of  the  TDOA  observations  is  low.  In  these  areas, 
the  ability  of  the  algorithm  to  estimate  accurately  the  location  of  the  emitter  is  limited. 
These  blind  areas  can  be  avoided  with  different  sensor  configurations  or  more  receivers. 

The  accuracy  of  the  algorithm  filter  is  also  heavily  dependent  on  choosing  the 
correct  values  for  the  covariance  of  the  state  excitation,  Q.  If  this  value  is  too  low,  the 
filter  reaches  steady  state,  but  the  steady  state  estimate  of  the  location  is  far  from  the  true 
location  of  the  emitter.  As  Q  is  increased,  the  algorithm  reaches  a  steady  state  around  the 
emitter's  true  location,  but  the  error  covariance  and  the  error  present  in  the  final  estimate 
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rises  much  higher.  Consequently  the  size  of  the  3a  error  ellipsoids  and  the  maximum 
error  present  in  the  final  estimates  are  much  larger. 

The  final  approaches  used  in  each  of  the  scenarios  yield  good  results.  In  the 
second  approach  used  in  scenario  #1,  the  accuracy  obtained  from  the  filtering  is 
approximately  3  degrees  in  bearing  and  about  10  percent  range.  For  scenario  #2  the  error 
is  approximately  5  degrees  in  bearing  and  5  percent  in  range,  and  in  the  third  scenario,  the 
error  is  approximately  2  degrees  in  bearing  and  3  percent  in  range. 

The  pulse  TDOA  algorithm  filters  the  TDOA  observations  one  at  a  time.  This 
filtering  technique  tends  to  skew  the  error  ellipsoids  and  align  them  with  the  loci  of 
constant  TDOA  that  corresponded  to  the  last  TDOA  observation.  This  places  the  major 
axis  of  the  error  ellipsoid  along  this  loci  and  tends  to  increase  the  error  in  that  direction. 

If  both  of  the  TDOA  observations  are  proceswd  simultaneously  as  the  filter  reaches 
steady  state,  the  error  covariance  and  the  error  ellipsoids  can  be  decreased. 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Both  of  the  Kalman  filtering  algorithms  considered  here  perform  well.  Given  the 
proper  parameters,  both  algorithms  accurately  estimate  the  location  of  the  emitter  to  a 
reasonable  degree  of  accuracy.  The  orthogonality  of  the  TDOA  observations  is  the 
largest  factor  that  affects  the  accuracy  of  the  final  estimates. 

Overall,  the  geometry  of  the  burst  TDOA  filtering  problem  yields  TDOA 
observations  with  better  orthogonality.  Even  when  the  distance  to  the  emitter  is  much 
larger  than  the  distance  between  the  receivers,  the  orthogonality  of  the  TDOA 
observations  are  good.  In  the  burst  TDOA  filtering  problem,  the  orthogonality  of  the 
observations  is  not  as  sensitive  to  the  location  of  the  emitter  in  relation  to  the  receivers. 
The  receiver  configuration  provides  good  orthogonality  for  a  wide  range  of  emitter 
locations.  To  obtain  good  orthogonality  in  the  pulse  TDOA  filtering  problem,  the 
receivers  require  wide  separation. 

The  major  disadvantage  of  the  burst  TDOA  filter  is  the  slow  rate  at  which 
observ'ations  are  obtained.  Observations  are  obtained  at  the  scan  rate  of  the  emitter,  and 
for  the  scenarios  examined,  this  means  that  only  one  observation  is  made  per  second.  For 
the  simulations  pursued,  the  total  run  time  is  25  seconds,  but  reasonably  accurate 
estimates  of  the  location  are  available  within  about  1 0  observations.  The  number  of 
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observations  required  to  reach  steady  state  could  be  reduced  further  with  better  a  priori 
estimates. 

The  most  attractive  feature  of  the  pulse  TDOA  filter  is  the  rate  at  which  observations 
are  obtained.  Observations  are  available  at  the  rate  of  the  PRF  of  the  emitter,  and  for  the 
scenarios  pursued  here,  the  pulse  TDOA  filter  could  obtain  2000  observations  for  every 
observation  obtained  by  the  burst  filter. 

It  may  be  possible  to  combine  the  two  filtering  algorithms  and  utilize  the  best 
features  of  both.  As  stated  previously,  with  only  two  receivers  the  location  of  the  emitter 
can  not  be  determined  uniquely  from  the  burst  TDOA  observations.  However,  if  the 
receivers  are  equipped  with  multiple  sensors,  as  in  the  pulse  TDOA  filtering  problem, 
bearing  estimates  could  be  obtained  and  the  location  of  the  emitter  along  the  burst  TDOA 
locus  could  be  uniquely  determined. 

B.  RECOMMENDATIONS 

The  analysis  and  development  conducted  here  is  limited  in  its  scope.  Further  testing 
and  evaluation  are  required  to  more  accurately  evaluate  the  algorithms  and  models 

presented.  Specifically,  the  following  are  recommended 

1 .  Further  evaluation  and  testing  of  the  burst  and  pulse  TDOA  algorithms  be  done  to 
assess  their  performance  for  a  wider  variety  of  receiver  and  emitter  configurations 
and  scenarios. 

2.  Perform  a  detailed  analysis  of  the  effect  of  a  priori  estimates  on  the  convergence 
and  accuracy  of  both  filters. 

3.  Explore  the  possibility  of  using  the  orthogonality  of  the  TDOA  observations  to 
pick  the  best  observations  to  filter  and  its  effect  on  the  performance  of  the  filter. 

4.  Obtain  more  detailed  models  of  the  error  present  in  the  burst  and  pulse  TDOA 
observations,  and  determine  their  effect  on  the  performance  of  the  filter. 

5.  Explore  the  possibility  of  combining  the  two  filtering  approaches  to  take 
advantage  of  the  best  features  of  both. 
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APPENDIX  A.  TOA  PROBABILITY  DENSITY 
A.  PULSE  TOA  PROBABILITY  DENSITY 

1 .  Derivation  of  Probability  Density 

The  pulses  received  and  detected  by  the  sensors  are  assumed  to  have  the 
envelope  shown  in  Figure  A-1 . 


Figure  A-1.  Pulse  envelope  and  sampled  pulse 


The  sampling  interval  is  assumed  much  smaller  than  the  pulse  width  so  the  sampled 
pulse  is  a  reasonable  representation  of  the  original  pulse.  The  higher  the  sampling 
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probability  density  function  of  the  pulse  TO  A,  each  sample  of  the  pulse  is  considered  a 
random  variables  as  shown  in  equation  (A-1). 

2(k)  =  Zo(k)  +  v(k)  (A- 1 ) 

Where  z(k)  =  the  random  variable  representing  the  amplitude  of  the  kth  sample, 

Zg(k)  =  the  deterministic  variable  that  is  the  true  amplitude  of  the  kth 
sample,  and 

v(k)  =  the  white  Gaussian  random  variable  that  represents  the  noise 
present  in  the  amplitude  of  the  kth  sample.  The  noise  has  zero 
mean  and  a  variance  of  V,  i.e..  E[v(k)^]  =  V. 

The  mean  and  variance  of  z(k)  are; 

Pz(k)  =  E[zo(k)  +  v(k)]  =  E(Zo(k)]  +  E[v(k)]  =  Zo(k)  (A-2) 

=  El  (z(k)  -  P2(k))(z(k)  -  pz(k))  ]  (A-3) 

Ojoc)  =  E[  (zo(k)  +  v(k)  -  Zo(k))(z„(k)  +  v(k)  -  Zo(k))  ]  =  E[v(k)2]  =  V  (A-4) 
Each  sample  is  uncorrelated. 

E[z(k)z(k  +  t)]  =  E[ZoOt)Zo(k  +  x)]  =  0  (A-5) 

Since  they  are  Gaussian  random  variables,  they  are  independent. 
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Equation  (A-6)  is  used  to  calculate  the  time  of  arrival  for  the  centroid  of  the  pulse. 


Y  I  t(k)2(k) 

TOA  =  §  =  !^=i5 -  (A-6) 

Where  TOA  =  the  time  of  arrival  for  the  pulse  centroid, 

X  =  numerator  of  the  pulse  centroid  TOA  equation, 

Y  =  denominator  of  the  pulse  centroid  TOA  equation, 
t(k)  =  the  time  at  which  the  kth  sample  of  the  pulse  is  taken, 
z(k)  =  the  amplitude  of  the  kth  sample  of  the  pulse,  and 
N  =  the  number  of  samples  in  the  pulse. 

The  time  base  from  which  the  time  of  the  samples  are  measured  is  arbitrary,  but  is 
common  for  all  TOA  measurements.  This  value  is  assumed  deterministic  and  is  therefore 
known  exactly  for  each  sample. 

The  numerator  and  denominator  of  equation  (A-6)  are  themselves  Gaussian  random 
variables.  As  shown  in  reference  [1],  if  W  is  a  linear  combination  of  independent 
Gaussian  distributed  random  variables. 


W  =  ai  X(l)  +  a2X(2)  +  a3X(3)  + . a„X(n)  (A-7) 

then  W  i.s  a  Gaussian  distributed  random  variable  with  a  mean  and  variance  given  by; 

HW  =  ai|iX(l)  +  a2tlX(2)  +  a3tlX(3)  + . anttxcn)  (A-8) 

I 

I  Ow  =  ajO^l) +  a2CT^2)  ®3tT^3)  +  •  -antJ^n)  (A-9) 


Based  upon  these  relations  the  mean  and  variance  of  the  numerator  and  denominator  of 
equation  (A-6)  are  given  by; 
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(A-10) 


\ix  =  I  t(k)zo(k) 
k»l 


o'x=It(k)'V  (A-ll) 

k«i 

^Y=2^Zo(k)  (A- 12) 


o^=  I  V  =  NV 

k=l 


(A- 13) 


Equation  (A-6),  the  estimate  of  the  TOA  of  the  pulse  centroid,  is  the  ratio  of  two 
Gaussian  distributed  random  variables.  The  joint  probability  density  function  for  two 
Gaussian  distributed  random  variables  x  and  y  is; 


exp 


K\y)  = 


2(l-pi) 


2jtOxOy(l  -P^) 


2\l/2 


Where  p,  =  the  mean  of  x  and  y, 

o,  Oy  =  the  standard  deviation  of  x  and  y, 
p  =  the  correlation  coefficient  which  is  given  by; 


(A-14) 


P  = 


cov[xy]  E[(x-  px)(y-  py) 

CTxOy  OxOy 


(A-15) 


For  equation  (A-6),  the  covariance  and  correlation  coefficient  are  calculated  from  the 
following: 


cov[X,  Y]  =  E[(S  t(k)z(k)  -  i  t(k)zo(k))(i  z(k)  -  S  Zo(k))l  (A- 1 6) 

k«l  k=l  k=l  k=l 


cov[X,  Y]  =  E[(|^  t(k)v(k))(^I  v(k))] 


(A- 17) 
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Since  the  noise  term  v(k)  is  considered  white  noise  with  zero  mean  and  variance  V,  the 
cross  terms  in  the  multiphcation  in  equation  (A- 17)  will  be  zero  and  the  resulting 
covariance  is; 

cov[X,  Y]  =  t(k)v(k)2  =  t(k)V  (A- 1 8) 

Thus  the  correlation  coefficient  is: 

VEt(k)  Et(k) 

n  = - — - =  — -  rA-lQi 

The  correlation  coefBcient  does  not  depend  upon  the  variance  of  the  noise  V  only  on  the 
number  of  samples  and  the  sampling  interval. 

From  reference  [2],  the  probability  density  function,  f^(z),  for  the  raho  of  two  random 

y 

variables,  Z  =  ^,  is  found  from  the  following  equation: 

fz(z)=  I  |y|f(zy,y)dy  (A-20) 

-«0 

Where  f(zy,y)  is  the  joint  probability  density  function  f(x,y)  with  the  variables  zy 
substituted  for  x. 
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The  probability  density  flmction  for  the  pulse  centroid  TOA  is: 

i  -  2.a.a,(.-pV^ - ^ 


(A-21) 


Where  z  =  the  pulse  centroid  TOA, 

y  =  the  denominator  of  the  pulse  centroid  TOA  equation , 

=  the  mean  of  the  numerator  and  denominator  of  the  centroid  TOA 
equation,  calculated  from  equations  (A- 10)  and  (A- 12), 
a/  =  the  variance  of  the  numerator  and  denominator  of  the  centroid 

TOA  equation,  calculated  from  equations  (A-1 1)  and  (A-13),  and 
p  =  the  correlation  coefficient  calculated  from  equation  (A- 15), 

A  closed  form  solution  for  this  integral  is  not  available,  so  a  MATLAB  program  is  used  to 
numerically  calculate  the  probability  density  flmction.  The  MATLAB  program  Puldist.m, 
in  Appendix  E,  calculates  the  probability  density  function  for  the  pulse  centroid  TOA  . 

2.  Effect  of  Signal  to  Noise  Ratio  on  Probability  Density 

The  plots  in  Figure  A-2  demonstrate  the  effect  of  the  peak  signal  to  noise  ratio  on 
the  probability  density  flmction  of  the  pulse  TOA.  Samples  of  white  Gaussian  noise  are 
added  to  the  1 .0  microsecond  pulse  shown  in  Figure  A-1 .  The  variances  of  the  noise  are 
chosen  to  give  peak  signal  to  noise  ratios  of  IS,  20,  25,  and  30  dB.  The  variance  of  the 
noise  is  calculated  from  the  peak  signal  to  noise  ratio  with  equation  (A-22). 


(A-22) 


Where 


V  =  the  variance  of  the  noise  required  to  give  the  specified  peak 
signal  to  noise  ratio  for  the  unit  pulse. 
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Figure  A-2.  Pulse  TOA  probability  density  for  peak 
signal  to  noise  ratios  of  1 5,  20,  25,  and  30  dB 


As  expected,  as  the  signal  to  noise  ratio  increases,  the  variance  of  the  error  present  in  the 
estimate  of  the  TOA  decreases  and  the  estimate  of  the  pulse  TOA  becomes  more  accurate. 

3.  Calculated  Mean  and  Standard  Deviation  of  Sampled  Pubes 

The  mean  and  standard  deviation  are  calculated  for  the  one  microsecond  pulse 
shown  in  Figure  A-1 .  A  variety  of  peak  signal  to  noise  ratios  and  sampling  rates  are  used. 
For  all  of  the  sampled  pulses  the  mean  time  of  arrival  for  the  centroid  is  0.5  microseconds. 
The  calculated  standard  deviations  are  shown  in  Table  A-I . 
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Sampling  Rate 

S/N 

4MHz 

8  MHz 

12MHz 

l6MHz 

15  dB 

005873 

0  02910 

0  02241 

0  01883 

20  dB 

0  03221 

0  01629 

001258 

001058 

25  dB 

001819 

0  00925 

000713 

000598 

30  dB 

0  01062 

0  00542 

0  00412 

000344 

TABLE  A-l.  STANDARD  DEVIATION  FOR  1  0  MICROSECOND 
SAMPLED  UInTIT  PULSE  (microseconds) 


The  pulse  TOA  standard  deviations  listed  in  Table  A-l  indicate  that  there  is  a 
tradeoff  between  sampling  rate  and  peak  signal  to  noise  ratio  If  a  lower  sampling  rate  is 
used,  a  larger  peak  signal  to  noise  ratio  is  required  to  obtain  the  same  TOA  error  level 
There  exists  an  mversely  proportional  relationship  between  the  variance  of  the  pulse  TOA 
error  and  the  peak  signal  to  noise  ratio  From  Table  A>l,  the  relationship  between  the 
variance  of  the  TOA  error  and  the  peak  signal  to  noise  ratio  for  a  sampling  rate  of  8  MHz 

Van^  =  (0.02910)^  ^  ^  _ 

Var2MB  (0.00925)2  '  ^ 


B.  BURST  TOA  PROBABILITY  DENSITY 

1.  Derivation  of  Probability  Density 

The  algorithm  that  calculates  the  probability  density  function  for  the  sampled 
pulse  is  used  to  calculate  the  probability  density  function  for  a  burst  of  pulses.  The 
MATLAB  program  Burdist.m,  listed  in  Appendix  E,  generates  the  burst  shown  in 
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Figure  A-3  and  calculates  the  p  ^ility  density  function  for  a  variety  of  peak  signal  to 
noise  ratios,  PRFs  and  pulse  wiuuis 


Figure  A-3.  A  burst  of  sampled  pulses  from  a  scanning  emitter 

The  envelope  of  the  burst  is  approximated  as  the  upper  half  of  a  sinewave,  and 
the  peak  signal  to  noise  ratios  are  calculated  using  equation  (A-22). 

2.  Effect  of  Signal  to  Noise  Ratio  on  Burst  TOA  Probability  Density 

Figure  A-4  demonstrates  the  effect  of  signal  to  noise  ratio  on  probability  density 
of  the  burst  TOA.  The  probability  densities  for  the  burst  in  Figure  A-3  are  calculated  for 
peak  signal  to  noise  ratios  of  1 5,  20,  2S,  and  30  dB.  The  noise  variances  required  are 
calculated  using  equation  (A-22). 


89 


Figure  A-4.  Burst  TOA  probability  densities  for  a  2  millsecond  burst 
with  a  pulsewidth  of  1 .0  miCToseconds,  and  a  PRF  of  6  kHz. 

As  expected,  as  the  signal  to  noise  ratio  is  increased,  the  variance  of  the  error 
present  in  the  estimate  of  the  TOA  decreases  and  the  estimate  of  the  burst  TOA  becomes 
more  accurate. 

3.  Effect  of  PRF  on  the  Burst  TOA  Probability  Density 

The  probability  density  is  calculated  for  a  burst  with  a  variety  of  pulse  repetition 
frequencies  (PRF).  The  PRFs  chosen  are  3, 6,  9,  and  12  kHz.  Plots  showing  the  bursts 
used  are  presented  in  Figure  A-S. 
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PRF  =  3  kHz  PRF  =  6  kHz 


PRF  =  9kHz 


PRF  =  12  kHz 


Figure  A-5.  Bursts  used  to  examine  the  efifect  of  PRF 
on  the  TOA  probability  density 


The  other  burst  characteristics  are: 

Burst  Length 
Pulsewidth 

Peak  Signal  to  Noise  Ratio 
Sampling  Rate 


2  0  milliseconds 
1 .00  microseconds, 
25  dB, 

8  MHz 
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The  calculated  probability  density  functions  for  the  bursts  in  Figure  A-S  are  plotted  in 
Figure  A-6. 


Figure  A-6.  Burst  TOA  probabibty  densities  for 
PRFs  of  3, 6,  9,  and  12  kHz 


As  the  PRF  increases  the  standard  deviation  of  the  burst  TOA  error  decreases. 
With  an  increase  in  the  PRF  more  energy  is  contained  in  the  burst  and  the  effect  of  the 
noise  present  diminishes 

4.  Calculated  Mean  and  Standard  Deviation  of  Sampled  Bursts 

The  means  and  standard  deviations  are  calculated  for  the  bursts  in  Figure  A-S  for 
a  variety  of  peak  signal  to  noise  ratios.  The  mean  time  of  arrival  for  the  centroid  of  the 
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burst  is  1 .0  milliseconds  for  all  of  the  bursts.  The  standard  deviations  are  Usted  in 
Table  A-2. 


Sampling  Rate 


S/N 

3000  Hz 

6000  Hz 

9000  Hz 

12000  Hz 

15dB 

0.02392 

0.01665 

0.01353 

0.01185 

20  dB 

0.01337 

0.00925 

0.00755 

0.00655 

25  dB 

0.00754 

0.00524 

0.00428 

0.00373 

30  dB 

0.00428 

0.00300 

0.00248 

0.00218 

TABLE  A-2.  STANDARD  DEVIATION  FOR  2.0  MILLISECOND 
SAMPLED  UNIT  BURST  (milliseconds) 


The  data  in  Table  A-2  show  that  the  PRFs  and  the  variances  of  the  burst  TOA  are 
inversely  proportional.  The  relationship  between  the  PRF  and  the  variance  of  the  burst 
TOA  is  shown  below  for  a  signal  to  noise  ratio  of  20  dB; 

Var,2ooo  ^  (0  00655)^  -  0  240  «  3000  Hz 
Varjooo  (0.01337)^"  12000  Hz 

Similar  to  the  relationship  between  the  PRF  and  the  variance  of  the  burst  TOA,  the  peak 
signal  to  noise  ratio  and  the  variance  of  the  burst  TOA  are  inversely  proportional.  For  a 
PRF  of  6000  Hz 


VarndB 

VarjsdB 


(0  01665)^ 
(000524)^ 


=  10  0s(25dB-  I5dB) 
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5.  Effect  of  the  Pulsewidth  on  the  Probability  Density 


The  probability  density  function  of  the  burst  TOA  is  calculated  to  er.amine  the 
effect  of  th;  pulsewidth  on  the  probability  density.  A  burst  is  generated  with  the  following 
characteristics; 

Burst  Length;  2.0  milliseconds 

Peak  Signal  to  Noise  Ratio  20  dB, 

Pulse  Repetition  Frequency  6000  Hz, 

Sampling  Rate  8  MHz 


For  this  burst  the  pulse  width  is  varied  from  1 .0  microseconds  to  4.0 
microseconds  and  the  probability  density  function  is  calculated  and  the  mean  and  standard 
deviations  are  examined.  The  probability  densities  for  each  of  the  bursts  are  plotted  in 
Figure  A-7. 
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0.95  1.0  1.05 


Centroid  Time  Of  Arrival  (milliseconds) 


Figure  A-7.  Burst  TOA  probability  densities  for 
pulse  widths  of  1.0,  2.0,  3.0,  and  4.0  microseconds 


As  shown  in  Figure  A-7,  the  effect  of  the  pulsewidth  on  the  probability  density 
function  is  similar  to  the  effect  of  the  PRF  and  the  peak  signal  to  noise  ratio  on  the 
probability  density  function.  As  the  pulsewidth  increases,  the  standard  deviation  of  the 
probability  density  decreases.  As  the  pulse  width,  and  effectively  the  power  present  in  the 
burst  increases,  the  estimate  of  the  TOA  of  the  centroid  of  the  burst  becomes  more 
accurate.  The  pulsewidth  and  the  variance  of  the  TOA  of  the  burst  centroid  are  inversely 
proportional.  The  relationship  between  the  pulse  width  and  the  variance  of  the  probability 
density  function  is  shown  below  for  a  signal  to  noise  ratio  of  20  dB  and  a  PRF  of  6000 
Hz: 
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Var4.o  _  (0.00472)^  _  ^  ^  1.0  ^isec 

Varj.o  (0.00926)^  4.0  ^  sec 
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APPENDIX  B.  LOCI  OF  CONSTANT  TDOA 


A.  BURST  TDOA  PROBLEM 

The  burst  TDOA  for  two  widely  separated  receivers  being  scanned  by  a  constantly 
rotating  emitter  is  directly  proportional  to  the  angle  formed  by  the  two  receivers  and  the 
emitter. 


y  Emitter 


Figure  B-1.  Angular  relationships,  burst  TDOA  problem 
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Using  the  geometry  in  Figure  B-1,  a  relationship  exists  between  the  bearing  and  range 


from  receiver  A  to  receiver  B,  the  bearing  and  range  to  the  emitter,  and  the  angle  6. 

_ Ra _ Rab 

sin(7i  -  (0  +  4)  -  \|/))  sin(6)  ^ ) 

Where  SR  =  the  scan  rate  of  the  emitter  in  radians/second, 

Ra  =  the  range  to  the  emitter, 

4*  =  the  bearing  to  the  emitter, 

Rb  =  the  range  to  receiver  B, 

\|/  =  the  bearing  to  receiver  B,  and 
6  =  the  angle  formed  by  the  receivers  and  the  emitter. 

When  this  equation  is  rearranged; 


tan(0)  = 


sin(4)-\|/) 

SR[^"Cos(4>-m/)] 


(B-2) 


If  the  assumption  is  made  that  the  angle  6  is  small,  then  the  tangent  of  6  is 
approximately  equal  to  the  angle  6  and  the  TDOA  can  be  found  from  the  following 
equation: 


TDO,.,-  e  _  sin((|i-v) 

SR  SR[ft-cos(*-v)] 


(B-3) 


Equation  (B-3)  is  used  to  plot  all  of  the  possible  locations  that  an  emitter  could  be  located 
given  a  TDOA  observation. 

B.  PULSE  TDOA  PROBLEM 

A  simple  relationship  exists  that  can  be  used  to  plot  the  loci  of  constant  TDOA  for  the 
pulse  TDOA  problem.  The  pulse  TDOA  problem  geometry  is  shown  in  Figure  B-2. 
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Emitter 


y 


Figure  B-2.  Pulse  TDOA  problem  geometry 

From  the  law  of  cosines; 

Rb^  =  Ra^+Rab^-2RaRabcos(<|)-\|/)  (B-4) 

The  difiTerence  between  Ra  and  Rb  is  the  distance  that  a  pulse  must  travel  to  reach  the 
more  distant  receiver  A.  The  amount  of  time  required  to  travel  this  distance  is  the  pulse 
TDOA.  This  distance  D  can  be  found  from  the  following  equation; 

D  =  (TDOA)c  (B-5) 

Where  TDOA  =  the  pulse  time  difference  of  arrival, 

c  =  the  speed  of  light. 


99 


Therefore,  Rb  is: 


Rb  =  Ra-D 


(B-6) 


Substituting  into  equation  (B>4): 

(Ra  -  D)2  =  Ra^  +  Rab^  -  2RaRab  cos(4  -  h;)  (B-7) 


<t)  =  arccos 


'Rab^-2RaD  +  D^ 
2RaRab 


(B-8) 


Equation  (B-8)  can  be  used  to  plot  the  possible  locations  of  an  emitter  given  the  locations 
of  the  sensors  and  the  TDOA. 


APPENDIX  C.  TDOA  ORTHOGONALITY 


A.  BURST  TDOA  ORTHOGONALITY 

Equation  (C- 1)  is  an  approximate  reiationslup  used  to  calculate  the  loci  of  constant 
TDOA  for  the  burst  TDOA  problem 


TDOA  = 


sin(^-V) 

SR[^-cos(4>-\(// 


(C-1) 


Where  SR  =  the  scan  rate  of  the  emitter  in  radians/second, 

Ra  =  the  range  to  the  emitter, 

((>  =  the  bearing  to  the  emitter, 

Rb  =  the  range  to  receiver  B, 

\|/  =  the  bearing  to  receiver  B,  and 
6  s  the  angle  formed  by  the  receivers  and  the  emitter. 

A  measure  of  the  orthogonality  of  the  TDOA  observations  can  be  found  by  taking  the 


dot  product  of  the  unit  vectors  tangent  to  the  loci  of  constant  TDOA  at  the  coordinates  of 


the  emitter.  These  vectors  can  be  calculated  from  the  slope  of  the  line  tangent  to  the  loci 


at  the  point  of  interest.  The  dot  product  of  the  unit  vectors  tangent  to  the  loci  will  yield 
the  cosine  of  the  angle  formed  by  the  two  loci.  The  equation  for  the  dot  product  is  given 
below: 


Rt  1  •  Rt2  =  cos  0  (C-2) 

Where  Rt,  =  The  unit  vector  tangent  to  the  loci  number  1, 

Rt,  =  the  unit  vector  tangent  to  loci  number  2,  and 
0  =  the  angle  formed  by  the  two  vectors. 

The  cosine  of  the  angle  formed  by  the  vectors  tangent  to  the  loci  is  an  indication  of  the 
orthogonallity  of  the  loci.  If  the  cosine  of  the  angle  formed  by  the  loci  is  close  to  one  then 
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the  loci  are  nearl>  parallel  If  the  cosine  of  6  is  zero  the  loci  are  orthogonal  and  lie  at  right 
an^es  to  each  other  The  slope,  dyt/dxt,  of  the  loci  of  constant  TDOA  was  found  by 
expressing  equation  (C- 1 )  in  Cartesian  coordinates  and  implicitly  diflferentiating  An 
assumption  is  made  to  simplify  the  problem.  If  the  distance  from  receiver  A  to  the  emitter 
is  much  larger  than  the  distance  to  receiver  B,  the  ratio  Ra/Rab  will  be  much  larger  than 
one,  and  equation  (C-I)  can  be  approximated  with  equation  (C-3). 


TDOA  = 


Rabsin((j)-v|/) 

SR(Ra) 


(C-3) 


For  a  given  loci  the  TDOA,  Rab,  and  SR  will  be  constant  and  are  not  a  function  of  the 
location  of  the  emitter.  The  equation  can  be  rewritten  in  the  following  form; 


SR(TDOA)  sinil) cos sin M/ cost}) 

Rab  ~  Ra  ^ 

The  following  substitutions  are  made  to  solve  for  the  equation  in  terms  of  the  location  of 

the  emitter  xt  and  yt. 

sin  4)  =  ^  cos<t>  =  ^  Ra  =  [(xt-xa)^ +(yt-ya)^]'^  (C-5) 

Ke  Rfl 

The  resulting  equation  is; 


SR(TDOA)  _  (yt  -  ya)cos  xa)sin  \\i 

Rab  [(xt  -  xa)^  +  (yt  -  ya)^] 

The  cosv|/  and  the  siny  terms  are  retained  to  make  the  derivations  more  clear.  Implicit 
differentiation  was  used  to  calculate  the  derivative  dyt/dxt  of  equation  (C-6).  The  term  on 
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the  left  side  of  equation  (C-6)  is  a  constant  and  therefore  its  derivative  is  zero.  The  final 
derivative  is  given  by  equation  (C-7). 


dxt  [(yt  -  ya)^cos  \\i  -  2(xt  -  xaKyt  -  ya)sin  vj/  -  (xt  -  xa)^cos  v|/] 


Where  xt,  yt  =  the  x  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  X  and  y  coordinates  of  receiver  A,  and 
V  =  the  bearing  fi-om  recover  A  to  receiver  B. 

The  bearing  fi’om  receiver  A  to  receiver  B  is  calculated  fi'om  the  following  equation; 


\^  =  arctan 


(C-8) 


The  unit  vector  tangent  to  the  loci  of  constant  TDOA  at  the  coordinates  of  the  emitter  is 
given  by  the  following  equation; 


Where  Rt  -  the  unit  vector  tangent  to  the  loci,  and 

m  -  the  slope  of  the  loci  dyt/dxt  at  the  emitter. 


(C-9) 


B.  PULSE  TDOA  ORTHOGONALITY 

The  dot  product  of  the  unit  vectors  tangent  to  the  loci  of  constant  TDOA  give  the 
cosine  of  the  angle  between  the  loci  of  constant  TDOA  and  a  good  indication  of  the 
orthogonality  of  TDOA  observations.  The  unit  vectors  tangent  to  the  loci  are  calculated 
from  the  slope  of  the  loci.  The  slope  of  the  loci  of  constant  TDOA  cai^  be  found  by  taking 
the  derivative  dyt/dxt  of  equation  (C-10). 
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";TX)A  = 


(xt  -  xa)^  +  (yt  -  ya)^  1  -  [(xt  -  xb)^  +  (yt  -  yb)^ 


(C-10) 


Where  xt,  yt  =  the  x  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  x  and  y  coordinates  of  sensor  A, 

xb,  yb  =  the  X  and  y  coordinates  of  sensor  B,  and 
c  =  the  speed  of  li^t. 

An  assumption  will  simplify  the  problem  and  yield  an  approximate,  but  simpler  solution. 

If  the  distance  from  the  sensors  to  the  emitter  is  much  larger  than  the  distance  between  the 
sensors,  the  approaching  pulses  can  be  considered  plane  waves,  and  the  TDOA  for  the 


pulses  arriving  at  the  sensors  can  be  calculated  from  the  geometry  in  Figure  C-1 . 
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The  TE>OA  calculated  from  the  geometry  in  Figure  C~  I . 


TDOA  = 


Rab  cos(^  -  v) 


(C-11) 


Where  Rab  =  the  distance  from  sensor  A  to  sensor  B, 

=  the  bearing  to  the  emitter, 

\|/  =  the  bearing  from  sensor  A  to  sensor  B,  and 
c  =  the  speed  of  Ught. 

Using  trigonometric  identities  and  the  substitutions  in  equations  (C-12)  and  (C-13), 
equation  (C-1 1)  can  be  rewritten  as  equation  (C'14). 


cos4>  = 


_ (xt-xa) _ 

[(xt  -  xa)^  +(yt-  ya)^] 


(C-12) 


sin4>  = 


(yt-ya) 

[(xt  -  xa)^  +  (yt  -  ya)^] 


(C-13) 


c(TDOA)  (xt  -  xa)cos  m/  (yt  -  ya)sin  vy 
"  ((xt-xa)2+(yt-ya)2]"^ 


(C-M) 


Implicit  differentiation  is  used  to  calculate  the  derivative  of  equation  (C-1 4).  The 
derivative  of  the  left  hand  side  is  zero  because  for  a  givoi  loci  the  terms  on  the  left  hand 
side  of  equation  (C-1 4)  are  constant.  The  equation  for  the  slope  of  the  loci  of  constant 
TDOA  is  given  in  equation  (C-1 5). 


dyt  _  (yt  -  ya)^cos  \|/  -  (xt  -  xa)(yt  -  ya)sin  \|/ 
dxt  (xt-xa)(yt-ya)cos>|>-(xt-xa)^sinM/ 

Where  xt,  yt  =  the  x  and  y  coordinates  of  the  emitter, 

xa,  ya  =  the  x  and  y  coordinates  of  sensor  A 

xb,  yb  =  the  X  and  y  coordinates  of  sensor  B,  and 
\|/  =  the  bearing  from  sensor  A  to  sensor  B. 


(C-1 5) 
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The  bearing  from  sensor  A  to  sensor  B  is  calculated  from  the  following: 


«|/  =  arctan 


r(yb-ya)' 

[(xb-xa). 


(C-16) 


The  unit  vector  tangent  to  the  loci  of  constant  TDOA  at  the  coordinates  of  the  emitter  is 
given  by  the  following  equation: 


(C-17) 


Where 


Rt  =  the  unit  vector  tangent  to  the  loci,  and 
m  =  the  slope  of  the  loci  at  the  emitter. 


APPENDIX  D.  BURST  H**  DERIVATION 


The  equation  that  calculates  the  TDOA  observations  as  a  function  of  the  receiver 
locations,  the  emitter  locations,  and  the  scan  rate  of  the  emitter,  is  given  by; 


TBOA  =  ^ 


_ (xt  -  xa)(yt  -  yb)  -  (yt  -  yaXxt  -  xb) _ 

[(xt  -  xa)^  +  (yt  -  ya)^l  '^[(xt  -  xb)^  +  (yt  -  yb)^] 


(D-1) 


Where  xt,yt  =  the  x  and  y  coordinate  of  the  emitter, 

xa,  ya  =  the  X  and  y  coordinate  of  receiver  A, 

xb,  yb  =  the  X  and  y  coordinate  of  receiver  B,  and 
SR  =  the  scan  rate  of  the  emitter  in  rad/sec. 

The  following  MATLAB  function  was  used  to  evaluate  the  partial  derivative  of  the 


observation  equation  (D-1).  The  terms  and  derivatives  in  the  MATLAB  function  were 


calculated  using  a  symbolic  processor  contained  in  the  software  package  MATHCAD. 


function  [Hx,Hy]  =  hk3(xt,yt,xa,ya,xb,yb,SR); 

%  A  function  to  calculate  the  derivative  of  the  observation  equation 
%  This  ftmction  assumes  that  all  of  the  receivers  are  allowed  to  move,  and 
%  Their  locations  are  known.  The  scan  rate  is  known  a  priori. 

% 

%  Richard  W.  Williamson 
%  Dale:  18  February  1994 
%  Date  Revised:  18  July  1994 
% 

%  This  program  evaluates  the  partial  derivatives  of  equation  (D-1 ) 


%  Define  extra  variables  for  clarity. 

xa2  =  xa''2; 
ya2  =  ya''2; 
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xb2  ®  xb''2; 
yb2  *  yb''2; 


xt2  =  xt^2; 
xt3  =  xt'^S; 

yt2  =  yt''2; 
yt3  =  yt'^S; 

%  Denominator  of  the  original  TDOA  fimction;  Equation  G^-1) 
G  =  ( (xt-xa)^2  +  (yt-ya)'^2  )*(  (xt-xb)''2  +  (yt-yb)''2  ); 


dgdx  =  -4*xt*yt*yb  +  8*xt*xa*xb  -  2*xa*xb2  -  2*xa*yt2  -  2*xa*yb2 
2*xa2*xb  -  2*yt2*xb  -  4*yt*ya*xt  -  2*ya2*xb  +  4*xt3  - ... 
6*xt2*xb  +  2*xt*xb2  +  4*yt2*xt  +  2*xt*yb2  -  6*xt2*xa  +... 
2*xa2*xt  +  2*ya2*xt  +  4*xa*yt*yb  +  4*yt*ya*xb; 


dgdy  =  -2*xt2*yb  -  4*xt*xa*yt  -  2*xa2*yb  -  4*yt*xt*xb  -  2*ya*xt2  - ... 
2*ya*xb2  +8*yt*ya*yb  -  2*ya*yb2  -  2*ya2*yb  +  4*yt3  +... 

4*xt2*yt  +  2*xa2*yt  +  2*yt*xb2  -  6*yt2*yb  +  2*yt*yb2  - ... 

6*yt2*ya  +  2*ya2*yt  +  4*xt*xa*yb  +  4*ya*xt*xb; 

%  Numerator  of  the  original  TDOA  fimction;  Equation  (D-1)  expanded  out 
f  =  xt*(ya-yb)  +  yt*(xb-xa)  +  (xa*yb-xb*ya); 

dfdx  =  (ya-yb); 
dfdy  =  (xb-xa); 

%  The  calculated  numerators  of  the  derivative  fimctions 

numdhdx  =  G*dfdx  -  0.5*f*dgdx; 
numdhdy  =  G*dfdy  -  0.5*f*dgdy; 


%  Evaluate  the  partial  derivatives 

Hx  =  (l/SR)*numdhdx/(sqrt(G^3)); 
Hy  =  (l/SR)*numdhdy/(sqrt(G^3)); 
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APPENDIX  E.  MATLAB  PROGRAMS 


%  Burst.m 
% 

%  Name:  Richard  W.  Williamson 
%  Date:  18  Feburary  1994 
%  Date  Revised:  18  July  1994 
% 

%  This  program  implements  an  extended  Kalman  filter  to  estimate 
%  the  location  of  an  emitter  from  the  TDOAs  of  a  burst  between  three 
%  receivers.  The  position  of  the  receivers  are  known  exactly. 

%  The  algorithm  generates  noise  free  TDOA  obsenrations  and  adds  white 
%  Gaussian  noise  to  the  observations.  The  variance  of  the  white  Gaussian 
%  noise  is  specified  in  the  program  and  can  be  different  for  each  receiver. 

%  Twenty-five  observations  are  processed. 

% 

%  The  sensors  are  initially  located  at  the  following  coordinates: 

% 

%  Receiver  A  (500, 100) 

%  Receiver  B  (20000,  -2000) 

%  Receiver  C  (-20000,-20000) 

% 

%  The  TDOA  measurements  are  corrupted  by  white  Gaussian  noise. 

%  The  variance  of  the  noise  for  all  TDOA  observations  is  1.694e-15. 

% 

%  The  emitter  was  located  at  a  range  of  500  kilometers  from  the  origin 
%  and  at  a  bearing  of  90  degrees  measured  from  the  x-axis. 

%  The  coordinates  of  the  emitter  were: 

% 

%  Emitter:  (0,  500)  kilometers. 

% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 
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% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


The  program  outputs  four  figures. 

figure(1)  The  plot  of  the  estimates  of  the  x  and  y  coordinates  of  the 
emitter 

figure(2)  The  Kalman  gains  calculated  for  each  observation 
figure(3)  An  X-Y  plot  of  the  estimates  of  the  coordinates  of  the  emitter. 
This  plot  also  includes  the  loci  of  constant  TDOA  and  the 
3o  error  ellipsoids  for  the  1st,  8th,  15th,  and  22nd  estimate 
of  the  location. 

figure(4)  A  closeup  of  the  steady  state  estimate  of  the  emitter  location. 
The  plot  is  centered  on  the  steady  state  estimate  of  the 
emitter  location  and  the  3a  error  ellipsoid  corresponding 
to  the  steady  state  estimate  is  plotted. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
The  program  calls  four  subroutines. 

tdoa3.m  Calculates  the  time  difference  of  arrival  of  a  burst 

between  two  receivers.  The  locations  of  the  receivers 
are  arbitrary.  The  TDOA  is  calculated  from  the  locations 
of  the  receivers,  the  location  of  the  emitter,  and  the  scan 
rate. 

hk3.m  Calculates  the  derivative  of  the  observation  equation  given 
the  locationof  the  two  receivers,  the  scan  rate,  and  the 
estimated  location  of  the  emitter. 

bloci.m  Calculates  the  loci  of  constant  TDOA  given  the  TDOA 
observation  and  the  locations  of  the  receivers. 

ellip.m  Calculates  the  error  ellipsoids  given  the  location  of  the  point 
of  interest,  the  error  covariance  matrix,  and  the  number  of 
standard  deviations  required. 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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clear; 


%  Range  and  bearing  to  the  target 
Rt  =  500000; 

Bmgt  =  90; 

xt  s  Rt*cos(Bmgt*pi/180); 
yt  =s  Rt*sin(Bmgt*pi/180); 

%  The  number  of  observations  processed. 

k=  1:25; 
len  =  length(k); 

%  The  coordinates  of  receivers  A,  B,  and  C  as  a  function  of  k 

xa  =  500*ones(1,len)  +  00*k;; 

ya  =  100*ones(1,len)  +  30*k;; 

xb  =  20000*ones(1.len)  +  30*k; 

yb  =  -2000*ones{1.len)  -  30*k; 

xc  =  -20000*ones(1,len)  -  30*k; 

yc  =  "20000*ones(1  ,len)  -  30*k; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

%  Initializing  variables 

Po=  1.00e+20*eye(2); 

Q  =  4.0e+06*eye(2); 

RKab=  5.660e-09; 

RKac=  5.660e-09; 

RKbc=  5.660e-09; 

Srate  =  2*pi; 

X  =  zeros(2,len); 

Pk  =  zeros{2,2); 

HK  =  zeros(3,2); 
g  =  zeros(2,len); 
z_tru  =  zeros(3,len); 

Z  =  zeros(3.len); 

Zhat  =  zeros(3,len); 


%  The  initial  error  covariance  matrix 
%  Variance  of  the  state  excitation  noise. 

%  The  variance  on  the  TDOA  error  Rcvr  A-B 
%  The  variance  on  the  TDOA  error  Rcvr  A-C 
%  The  variance  on  the  TDOA  error  Rcvr  A-C 
%  The  scan  rate  of  the  emitter 

%  The  state  i  iatrix 
%  The  error  covariance  matrix 
%  The  linearized  derivative  of  h(k) 

%  The  Kalman  gains 
%  The  vector  of  true  observations 
%  The  vector  of  noisy  observations 
%  The  vector  of  estimated  obsen/ations 
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%  Calculation  of  the  observations  noise  free  observations 
fork*  1:len, 

z_tru(1,k)  =  tdoa3(xt,yt,xa(k),ya(k),xb(k),yb(k),Srate); 
2_tru(2,k)  =  tdoa3(xt,yt,xa(k).ya(k),xc{k),yc(k),Srate); 
z_tru(3,k)  =  tdoa3(xt,yt,xb(k),yb(k),xc(k),yc(k).Srate); 
end; 

%  Calculate  the  noise  on  the  TDOA  estimates 

rand(’normar): 

n  =  sqrt(RKab)*rand(3,len); 

%  Form  the  matrix  of  noisy  observations. 

Z  =  2_tru  +  n; 

%  A  priori  estimate  of  the  location  of  the  emitter 
Xo  =  [0  10000]': 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Initialize  the  Extended  Kalman  filter 

Pk  =  Po:  %  Initialize  the  Pk  the  error  covariance  matrix. 

X(:.1)  =  Xo;  %  Initialize  the  state  vector 

for  k  «  1:len, 

k 

%  Processes  the  observation  for  receiver  A  and  B 

^0  Calculates  the  derivative  of  the  Rcvr  A-B  TDOA  equation  wrt  x  and  y 
[hlab  h2ab]  =  hk3(X(1.k).X(2.k).xa(k).ya(k).xb(k),yb(k),Srate): 

%  Forms  the  Hk  matrix 
HKab  =  [hlab  h2ab]; 

%  Calculates  the  next  PK  based  upon  the  previous  PK 


Pk  =  Pk  +  Q: 
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%  Calculates  the  kalman  gains 

GKab  =  Pk*HKab'*jnv(  HKab*Pk*HKab’  +  RKab  ); 

%  Calculates  the  estimate  of  the  observation  based  upon  the 
%  estimate  of  the  states. 

Zhat(1,k)  =  tdoa3(X(1.k).X(2.k),xa(k).ya(k).xb(k).yb(k).Srate): 

%  Calculates  an  update  of  the  state  based  upon  the  difference  in 
%  the  actual  observation  and  the  estimate  of  the  observation. 
X(:,k)  =  X(:,k)  +  GKab*(  Z(1.k)-Zhat(1.k) ); 

%  Updates  the  variance  matrix  Pk  based  upon  the  observation 
Pk  =  (eye(2)  -  GKab*HKab)*Pk; 

%  Processes  the  observation  for  receiver  A  and  C 

%  Calculates  the  derivative  of  the  TDOA  equation  wrt  x  and  y 
[h1ac  h2ac]  hk3(X(1.k).X(2.k).xa(k).ya(k).xc(k).yc(k),Srate): 

%  Forms  the  Hk  matrix 
HKac  =  [hlac  h2ac]: 

%  Calculates  the  next  PK  based  upon  the  previous  PK 
Pk  =  Pk  +  Q: 

%  Calculates  the  kalman  gains 

GKac  =  Pk*HKac'*inv(  HKac*Pk*HKac’  +  RKac ); 

%  Calculates  the  estimate  of  the  observation  based  upon  the 
%  estimate  of  the  states. 

Zhat(2.k)  =  tdoa3(X(1,k),X(2,k),xa(k).ya(k).xc(k).yc(k).Srate); 

%  Calculates  an  update  of  the  state  based  upon  the  difference  in 
%  the  actual  observation  and  the  estimate  of  the  observation. 
X(:.k)  =  X(;,k)  +  GKac*(  Z(2.k)-Zhat(2.k) ); 

%  Updates  the  variance  matrix  Pk  based  upon  the  observation 
Pk  =  (eye(2)  -  GKac*HKac)*Pk: 
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%  Processes  the  observation  for  receiver  B  and  C 

%  Calculates  the  derivative  of  the  TDOA  equation  wrt  x  and  y 
[h1bc  h2bc]  =  hk3(X(1,k),X(2,k),xb(k).yb(k).xc(k).yc(k),Srate): 

%  Forms  the  Hk  matrix 
HKbc  =  [h1bc  h2bc]: 

%  Calculates  the  next  PK  based  upon  the  previous  PK 
Pk  =  Pk  +  Q: 


%  Calculates  the  kalman  gains 

GKbc  =  Pk*HKbc“inv(  HKbc*Pk*HKbc*  +  RKbc ); 

%  Calculates  the  estimate  of  the  observation  based  upon  the 
%  estimate  of  the  states. 

Zhat(3,k)  =  tdoa3(X(1,k),X(2,k).xb(k),yb(k),xc(k),yc(k),Srate): 

%  Calculates  an  update  of  the  state  based  upon  the  difference  in 
%  the  actual  observation  and  the  estimate  of  the  observation. 
X(;.k)  =  X(:.k)  +  GKbc*(  Z(3,k)-Zhat(3.k) ); 

%  Updates  the  variance  matrix  Pk  based  upon  the  observation 
Pk  =  (eye(2)  -  GKbc*HKbc)*Pk: 

%  Projects  the  next  state  based  upon  the  current  state 

X(;.k+1)  =  X(:.k): 

%  Save  the  error  covariance  matricies 
ifk==  1 
P  =  Pk; 
else, 

P  =  IPiPk]: 
end; 

%  Save  the  Kalman  gains 
g(;.k)  =  GKbc; 


end; 


%  End  of  the  filtering  routine 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%o^%o^% 
%  Output  the  filtered  estimates  of  x  and  y 


t  =  1  :len; 
figure(1) 

plot(t.X(1  .t)/1 000.t.X(2.t)/1 000) 
axis([1  25  -50  600]) 
pause; 

%  Plot  of  the  Kalman  gains 

figure(2) 

cig 

plot(t.g(1.:).t.g(2.:)) 
title('Kalman  Gains') 
xlabeK'Number  of  Observations  ) 
ylabel('Magnitude') 
pause 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Calculate  the  steady  state  estimate  of  emitter  location 

xtavg  =  mean(X(1,k-10:k)); 
ytavg  »  mean(X(2.k-10:k)); 


%  Calculate  the  coordinates  of  the  loci  of  constant  TDOA  based 
%  upon  noise  free  TDOA  observations 

Tab  =  bloci(xa(k),ya(k).xb(k).yb(k).z_tru(1.k),Srate); 

Tac  =  bloci(xa(k),ya(k).xc(k),yc(k),z_tru(2.k),Srate); 

Tbc  =  bloci(xb(k),yb(k),xc(k),yc(k),z_tru(3,k),Srate); 


%  Plot  the  locus  of  emitter  locations  and  actual  emitter  locations 

figure(3) 

cig 

%  Plots  coordinates  of  the  estimates  of  the  emitter  location 
%  axis([-500  500  -30  1000]); 
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plot(X(1  ,:)/1 000,X(2,:)/1  OOO.r ) 
hold  on 


%  Plots  the  actual  emitter  location 
plot(xt/1000.yt/1000.'*g') 

%  Plots  the  final  locations  of  the  receivers 
plot(xa(k)/1 000.ya(k)/1 000.’*g*) 
plot(xb(k)/1 000,yb(k)/1 000, '*g') 
plot(xc(k)/1 000.yc(k)/1 000.’*g') 

%  Plots  the  locus  of  possible  emitter  locations 
plot(  Tab(1.:)/1000,  Tab(2.:)/1000.’:b*) 
plot(  Tac(1,:)/1000.  Tac(2.:)/1000.':b’) 
plot(  Tbc(1.:)/1000.  Tbc(2.:)/1000.‘:b') 

%  Plot  three  sigma  error  elipsoids  for  various  points 

C  =  3.00;  %  The  number  of  standard  deviations  plotted  in  error  ellipsoids 

fork=  1:7:len"1, 

P1  *  [P(2*k-1 ,1 )  P(2*k-1 ,2);  P(2*k,1)  P(2*k,2)l:  %  The  Pk  Matrix 

lxg,y1  ,y2,Emax,Emin]  =  ellip(P1  ,C,X(1  ,k),X(2,k)):  %  Calculate  ellipse 

plot((xg)/1 000, (y1  )/1 000,'g-’,(xg)/1 000.(y2)/1 000,*g-’)  %  Plot  the  ellipse 

end; 

%  axis; 
hold  off 

axis([-250  250  -50  600]); 
pause 

%  Plot  the  locus  of  emitter  locations  and  actual  emitter  locations 

figure(4) 

cig 

%  Plots  coordinates  of  the  estimates  of  the  emitter  location 

plot(X(1  ,:)/1000,X(2,:)/1000,'r’ ) 
hold  on 
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%  Plots  the  actual  emitter  location 
plot(xt/1000.yt/1000;*g'): 

%  Plots  the  locus  of  possible  emitter  locations 
plot(  Tab(1.:)/1000.  Tab(2.:)/1000.‘:b'): 
plot(  Tac(1.:)/1000,  Tac(2.:)/1000.':b'); 
plot(  Tbc(1.:)/1000,  Tbc(2.:)/1000;:b'): 

%  Plots  the  error  ellipsoids 

[xg,y1,y2,Emax.Emin]  =  ellip(Pk.C.xtavg.ytavg); 

plot((xg)/1000.(y1)/1000;g-'.(xg)/1000.(y2)/1000;g-') 

Major  =  sprintf('%5.1f,Emax/1000): 

Minor  =  sprintf('%5.1f, Emin/1 000); 

text((xt-0.50*Emax)/1 000,(yt+0.85*Emax)/1 000, fMajor  Axis  =  ’  Major '  km*]) 
text((xt-0.50*Emax)/1 000, (yt+0.75*Emax)/1000, [Minor  Axis  =  *  Minor  ’  km*]) 

%  axis; 
hold  off 

axis([(xt-Emax)/1000  (xt+Emax)/1000  (yt-Emax)/1000  (yt+Emax)/1000]); 

%  axis('square') 

% 

% 

% 
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%  tdoa3.m 


% 

%  Name;  Richard  W.  Williamson 
%  Date:  18  Feburary  1994 
%  Date  Revised;  18  July  1994 
% 

%  This  program  calculates  the  burst  observation  equation  h(x(k)). 


% 

%  Input;  The  coordinates  of  receiver  A  (meters);  xa,  ya 

%  The  coordinates  of  receiver  B  (meters);  xb,  yb 

%  ITre  coordinates  of  the  emitter  (meters);  xt,  ^ 

%  Scan  rate  of  the  emitter:  SR 

% 

%  Output  The  burst  time  difference  of  arrival  for 
%  the  specified  receiver  loactions.  T 

% 

% 

function  T  =  TDOA3(xt,yt,xa,ya,xb,yb,SR) 


%  Numerator  of  the  observation  equation 
f  =  xt*(ya-yb)  +  yt*(xb-xa)  +  (xa*yb-xb*ya); 

%  Denominator  of  the  observation  equation 
G  =  ( (xt-xa)''2  +  (yt-ya)^2  )*(  (xt-xb)''2  +  (yt-yb)'^2 ); 

%  Calculate  the  TDOA 
T  =  f/(SR*sqrt(G)); 

% 
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%  HK3.m 
% 

%  Name:  Richard  W.  Williamson 
%  Date:  18  Feburary  1994 
%  Date  Revised:  18  July  1994 
% 


%  This  program  calculates  the  derivative  of  the  observation  equation  h(x(k)). 
% 


%  Input: 
% 

% 

% 

% 

%  Output 
% 

% 


The  coordinates  of  receiver  A  (meters):  xa,  ya 

The  coordinates  of  receiver  6  (meters):  xb,  yb 

The  coordinates  of  the  emitter  (meters):  xt,  ^ 

Scan  rate  of  the  emitter:  SR 

The  vector  of  partial  derivaties  of  h(x(k)) 

with  respect  to  xt  and  yt.  [Hx.Hy] 


function  [Hx.Hy]  =  hk3(xt,yt,xa.ya,xb,yb,SR); 


%  A  function  to  calculate  the  derivative  of  the 
%  H  matrix.  This  function  does  not  assume  that  one  of 
%  The  receivers  is  fixed  with  respect  to  the  stationary 
%  coordinate  system.  This  assumes  that  the  scan  rate  is  known. 
% 

%  This  program  is  an  implementation  of  equation  (  ) 

% 

%  Define  variable  for  easier  understanding 


xa2  =  xa''2; 
ya2  =  ya''2; 


xb2  =  xb''2: 
yb2  =  yb'‘2: 

xt2  =  xt''2: 
xt3  =  xt''3; 

yt2  =  yt''2; 
yt3  =  yt''3; 
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% 

%  The  original  TDOA  function: 


G  =  ( (xt-xa)''2  +  (yt-ya)''2  )*(  (xt-xb)^2  +  (yt-yb)''2 ); 

dgdx  *  -4*xt*yt*yb  +  8*xt*xa*xb  -  2*xa*xb2  -  2*xa*yt2  -  2*xa*yb2 
2*xa2*xb  -  2*yt2*xb  -  4*yt*ya*xt  -  2*ya2*xb  +  4*xt3  - ... 
6*xt2*xb  +  2*xt*xb2  +  4*yt2*xt  ^  2*xt*yb2  -  6*xt2*xa  +... 
2*xa2*xt  +  2*ya2*xt  +  4*xa*yt*yb  +  4*yt*ya*xb: 

dgdy  =  -2*xt2*yb  -  4*xt*xa*yt  -  2*xa2*yb  -  4*yt*xt*xb  -  2*ya*xt2  - 
2*ya*xb2  +8*yt*ya*yb  -  2*ya*yb2  -  2*ya2*yb  +  4*yt3  +... 
4*xt2*yt  +  2*xa2*yt  +  2*yt*xb2  -  6*yt2*yb  +  2*yt*yb2  - ... 
6*yt2*ya  +  2*ya2*yt  +  4*xt*xa*yb  +  4*ya*xt*xb; 

f  =  xt*(ya-yb)  +  yt*(xb-xa)  +  (xa*yb-xb*ya); 


dfdx  =  (ya-yb): 
dfdy  =  (xb-xa); 

%  The  calculated  numerators  of  the  derivative  functions 

numdhdx  =  G*dfdx  -  O.STdgdx; 
numdhdy  =  G*dfdy  -  O.STdgdy; 

%  The  actual  derivative  functions 

Hx  =  (1/SR)*numdhdx/(sqrt(G''3)); 

Hy  =  (1/SR)*numdhdy/(sqrt(G''3)): 
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Bloci.m 


% 


%  Name:  Richard  W.  Williamson 
%  Date:  18  Feburary  1994 
%  Date  Revised:  18  July  1994 
% 


% 

% 

% 

% 


% 

% 

% 

% 

% 


This  program  calculates  the  coordinates  of  the  loci  of  constant  TDOA  for 
a  pair  of  receivers. 


Input:  The  coordinates  of  receiver  A  (meters):  xa,  ya 

The  coordinates  of  receiver  B  (meters):  xb,  yb 

The  Time  Difference  of  Arrival  (sec):  TD 

Scan  rate  of  the  emitter:  SR 

Ou^ut  The  vector  of  x  and  y  coordinates  of  the  loci 

of  constant  TDOA.  T 


% 


% 

function  T  =  bioci(xa,va,xb,yb,TD,Srate) 


%  Polar  plot  of  the  receivers,  the  emitter, and  the  locus  of  possible 
%  emitter  location. 

%  The  angles  of  phi  to  be  plotted, 
phi  >  0:pi/100:pi; 


%  Calculates  the  range  and  bearing  to  the  receivers. 
Rab  =  sqrt(  (xb-xa)''2  +  (yb-ya)''2  ); 

psi  =  atan2(yb-ya,xb-xa)*ones(1,1:length(phi)); 

%  initialize  the  range  variables 
Ra  =  zeros(1,length(phi)); 

Ra  =  Rab*(  sin(phi-psi)Aan(TD*Srate)  +  cos(phi-psi) ); 

dxa  =  xa*ones(1,1:length(phi)); 
dya  =  ya*ones(1,1:length(phi)): 

%  Plots  the  locus  of  possible  emitter  location. 

T  =  [Ra.*cos(phi)+dxa  ;  Ra.*sin(phi)+dya  J 
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%  pulse,  m 
% 

%  This  program  implements  an  extended  Kalman  filter  to  estimate 
%  the  location  of  an  emitter  from  the  TDOA  of  a  pulse  between  two 
%  sensors.  Two  receiver  platforms  are  used  and  the  their  location  is 
%  known  exactly.  Each  receiver  platform  consists  of  two  sensors  The 
%  algorithm  generates  noise  free  TDOA  observations  for  the  known  receiver 
%  and  emitter  locations  and  adds  white  Gaussian  noise  to  the  observations. 
%  The  variance  of  the  white  Gaussian  noise  is  specified  in  the  program 
%  and  can  be  different  for  each  receiver  platform.  Sixty  obsevations 
%  are  processed. 

% 

% 

% 

%  The  sensors  are  located  at  the  following  coordinates: 

% 

%  Receiver  1  Sensor  A  (0,0) 

%  Receiver  1  Sensor  B  (500,0) 

%  Receiver  2  Sensor  A  (20000,0) 

%  Receiver  2  Sensor  B  (20500,0) 

% 

%  The  TDOA  measurements  are  corrupted  by  white  Gaussian  noise. 

%  The  variance  of  the  noise  for  both  receivers  is  1 .694e>1 5. 

% 

%  The  emitter  was  located  at  a  range  of  30  kilometers  from  the  origin 
%  and  at  a  bearing  of  60  degrees  measured  from  the  x-axis. 

%  The  coordinates  of  the  emitter  were: 

% 

%  Emitter:  (15,  26)  kilometers. 

% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 


% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 


The  program  outputs  four  figures. 

figure(1 )  The  plot  of  the  estimates  of  the  x  and  y  coordinates  of 
the  emitter 

figure(2)  The  Kalman  gains  calculated  for  each  observation 

figure(3)  An  X-Y  plot  of  the  estimates  of  the  coordinates  of 

the  emitter. 

This  plot  also  includes  the  loci  of  constant  TDOA  and  the 
3a  error  ellipsoids  for  the  first,  20th,  39th,  and  58th 
estimate  of  the  emitter. 

figure(4)  A  closeup  of  the  steady  state  estimate  of  the  emitter 
location.  The  plot  is  centered  on  the  steady  state 
estimate  of  the  emitter  location  and  the  3a  error 
ellipsoid  corresponding  to  the  steady  state  estimate 
is  plotted. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


The  program  calls  two  subroutines. 

ploci.m  Calculates  the  loci  of  constant  TDOA  given  the  TDOA 
observation  and  the  locations  of  the  receivers, 
ellip.m  Calculates  the  error  ellipsoids  given  the  location  of  the  point 
of  interest,  the  error  covariance  matrix,  and  the  number  of 
standard  deviations  required. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Initialize  the  variables 


c  =  3.00e+08; 

R1  =  1.694e-15: 

R2  =  1.694e-15; 

Q  =  5.0e+02*eye(2): 
PO  =  5.0e+20*eye(2): 


%  The  speed  of  light 

%  The  covariance  of  measurement  error  receiver  1 
%  The  covariance  of  measurement  error  receiver  2 
%  Variance  of  Plant  excitation  noise 
%  Initial  error  covariance 


%  Choose  to  calculate  the  error  covariance  of  3  sigma 
C  =  3.0, 


t  =  (0:60);  %  Time  vector  total  of  61  observations 

I  =  ones(1  ,len);  %  occurring  50  micro  seconds  apart 
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%  Initialize  storage  for  Kalman  gains  and  the  states, 
g  =  zeros(2,length(t)); 

X  =  zeros(2,len): 

%  Emitter  location 

Rt  -  30000;  %  Range  to  the  emitter 

Bt  =  60.00;  %  Bearing  to  the  emitter  from  the  origin 

xt  =  Rt*cos(Bt*pi/1 80);  %  X  coordinate  of  the  emitter 

yt  =  Rt*sin(Bt*pi/180);  %  Y  coordinte  of  the  emitter 

%  A  priori  estimate  of  the  location  of  the  target; 

XO  =  (10000;50001; 

% 

%  The  location  of  the  emitter  specified  for  all  point  in  time 

xa1  =  0*l  +  0^;  %  Receiver  1  Sensor  A  moves  in  a  straight  line 

ya1  =  0*l  +  0*t;  %  at  a  constant  speed 

xb1  =  500*l  +  0*t;  %  Receiver  1  Sensor  B  moves  in  a  straight  line 

yb1  =  0*l  +  0*t;  %  at  a  constant  speed 

xa2  =  20000*1  +  0*t;  %  Receiver  2  Sensor  A  moves  in  a  straight  line 

ya2  =  000.00*1  +  0*t;  %  at  a  constant  speed 

xb2  =  20500*1  +  0*t;  %  Receiver  2  Sensor  B  moves  in  a  straight  line 

yb2  =  000*1  +  0*t;  %  at  a  constant  speed 


%  Calculate  the  TDOA  observations  for  the  two  receivers 
%  Equation 

%  z(k+llk) 

%  Calculate  the  observations  for  receiver  1 

Z1  =  (1/c)*(  sqrt(  (xt-xa1).''2  +  (yt-ya1).''2  )  - ... 
sqrt(  (xt-xb1).''2  +  (>1-yb1).''2  ) ); 


%  Calculate  the  observations  for  receiver  2 
Z2  =  (t/c)*(  sqrt(  (xt-xa2).''2  +  (yt-ya2).''2  )  - ... 


(^)  -  xa(k))^  +  (>t(k)  -  ya(k))^]'^  -  \  (xt(k)  -  xb(k))^  +  (yt(k)  -  yb(k))^ 
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sqrt(  (xt-xb2).''2  +  (yt-yb2).^2  ) ); 


%  Inject  zero  mean,  noise  into  the  observations 

rand('normar): 

n1  =  sqrt(R1)*rand(1,len); 

n2  =  sqrt(R2)*rand(1,len); 

Z1n  =  Z1  +n1: 
iZ2n  =  Z2  +  n2: 

%  initialize  the  Kalman  filter 

X(; ,  1 )  :=  XO;  %  A  priori  estimate  of  the  emitter  location 

PK  =  PO;  %  Initial  error  covariance  matrix 

%  Calculate  the  filtered  estimates  of  the  emitter  location 

fork  =  1:len-1. 

%  Process  TDOA  observation  Receiver  #1 

/O  /O  A) /O  /OyO/O/D/OAl/O/DAt/vyO/O/O  /« /O  /D  /O  Al  Ai  /O  A>  VP /w  A>A>AlA>A>A>AiA>A) 


%  Calculate  the  HK  matrix  for  receiver  platform  1 

%  Equation  H/Tf'i  -I  U  _  (^)-x»>)"\  if  _  (^)-yb)"\ 

»  Ra(k)  Rb(k)  J  'V  Ra(k)  Rb(k)  J 

%  Equation  Ra(k)  =  [(xt(k)-xa)^ +  (yt(k)-ya)^]^^ 

%  Equation  Rb(k)  =  [(Jrt(k)  xb)^  +  (yt(k)  -  yby  ] 

Ra1  =  sqrt(  (X(1,k)-xa1(k)).''2  +  (X(2,k)-ya1(k)).''2  ); 

Rb1  =  sqrt(  (X(1,k)-xb1{k)).''2  +  (X(2.k)-yb1{k)).''2  ); 

h1x  =  (1/cr(  (X(1,k)-xa1(k))./Ra1  -  (X(1,k)-xb1(k))./Rb1  ); 
h1y  =  (1/c)*(  (X(2,k)-ya1(k))./Ra1  -  (X(2,k)-yb1(k))./Rb1  ); 


HK  =  [h1xh1y]; 


%  Calculate  the  estimates  of  the  TDOA  based  upon  the  estimate 
%  of  emitter  location  and  location  of  receiver  platform  1 
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%  Equation 
%  z(k+llk) 

Zhat1  =  (1/c)*(sqrt(  (X(1,k)-xa1(k))''2  +  (X(2.k)-ya1(k))''2  )  - ... 
sqrt(  (X(1.k)-xb1(k)r2  +  (X(2,k).yb1(k)r2  )); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q. 

PK  =  PK  Q; 

%  Calculate  the  Kalman  gains 
GK  =  PK*HK'*inv(HK*PK’‘HK  +  R1): 

%  Calculate  the  updated  error  covariance  matrix 
PK  =  (eye(2)  -  GK*HK)*PK; 

%  Calculate  a  smoothed  estimate  of  the  bearing  to  the  emitter 
X(;,k)  =  X(:,k)  +  GK*(Z1n(k)-Zhat1): 

%  estimate  of  the  emitter  location  for  receiver  #2  is  same  as  previous 

X(;,k)  =  X(;.k): 

%  Process  TDOA  observation  receiver  #2 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Calculate  the  HK  matrix  for  receiver  platform  2 

„/  ,,  ..  TT/i-N  r  iTcxtdchxa)  (3h(k)-xb)^  1  (^(k)-yb)'\ 

%  Equation  H(k)  =  - - 

%  Equation  Ra(k)  =  [(xt(k)  -  xa)^  +  (^(k)  -  ya)^  J  ^ 

%  Equation  Rb(k)  =  [(xt(k)  -  xb)^  +  (^(k)  -  yb)^  ]  ^ 

Ra2  =  sqrt(  (X(1,k)-xa2(k)).''2  +  (X(2.k)-ya2(k)).''2  ); 

Rb2  =  sqrt(  (X(1,k)-xb2(k)).^2  +  (X{2,k)-yb2(k)).''2  ); 
h2x  =  (1/c)*(  (X(1,k)-xa2(k))./Ra2  -  (X(1,k)-xb2(k))./Rb2  ); 
h2y  =  (1/c)*(  (X(2,k)-ya2(k))./Ra2  -  (X(2.k)-yb2(k))./Rb2  ); 


(^)  -  xa(k))^  +  (^)  -  ya(k))*  1  -  f  (xt(k)  -  xb(k))^  +  (^)  -  yb(k))^ 
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HK  =  [h2x  h2y]; 


%  Calculate  the  estimates  of  the  TDOA  based  upon  the  estimate 
%  of  the  emitter  location  and  location  of  receiver  Platform  2 


% 


% 


Equation 

z(k  +  1  Ik)  = 


(^)  -  xa(k))^  +  (yt(k)  -  ya(k))^ 


-  (^)  -  xb(k))^  4-  (^)  -  yfa(k))^ 


Zhat2  =  (1/c)*(sqrt(  (X(1.k)-xa2(k))''2  +  (X{2.k)-ya2(k))'^2  )  - ... 
sqrt(  (X(1.k)-xb2(k))'^2  +  (X(2.k)-yb2(k))>'2  )); 

%  The  new  error  covariance  is  the  same  as  the  old  plus  Q. 
PK  =  PK  +  Q; 


%  Calculate  the  Kalman  gains 
GK  =  PK*HK*inv(HK*PK*HK’  +  R2); 

%  Calculate  the  updated  error  covariance  matrix 

PK  =  (eye(2)  -  GK*HK)*PK: 

%  Calculate  a  smoothed  estimate  of  the  bearing  to  the  emitter 
X(:.k)  =  X(:,k)  +  GK*(Z2n(k)-Zhat2): 

%  Estimate  of  the  location  for  next  observations  is  same  as  previous 
X(;,k+1)  =  X(:,k); 

%  Save  the  error  covariance  matricies 
if  k==  1 
P  =  PK; 
else, 

P  =  [P;PK]; 
end; 

%  Save  the  Kalman  gains 
g(:.k)  =  GK; 

end;  %  end  of  the  Kalman  filtering  routine 
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%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

% 

%  Plot  the  output 
% 

% 


%  Plot  of  the  estimate  of  the  x  and  y  coordinates 

figure(1) 

cig; 

plot(t,X(1  ,;)/1 000.t.X(2.:)/1000) 
titleC  Estimate  of  x  and  y  coordinates'); 
xlabei('Number  of  Observations') 
ylabel('X  location  in  kilometers') 
pause 

%  Plot  of  the  Kalman  gains 

figure(2) 

cig 

plot(t.g(1.:).t.g(2.:)) 
title('Kalman  Gains') 
xlabeK'Number  of  Observations') 
yiabel('Magnitude') 
pause 

%  Calculate  the  steady  state  estimate  of  emitter  location 
xtavg  =  mean(X(1,k-10:k)); 
ytavg  =  mean(X(2,k-10:k)); 

%  Calculate  the  loci  of  constant  TDOA  based  upon  noise  free  TDOA 
%  observations 

R1a  =  5000:5000:40000; 

%  T1  is  a  vector  of  the  X  and  Y  coordinates  of  the  loci. 

T1  =  ploci(xa1(k),ya1(k),xb1(k),yb1(k).Z1(k).R1a); 

R2a  =  5000:5000:40000; 

%  T2  is  a  vector  of  the  X  and  Y  coordinates  of  the  loci. 

T2  =  ploci(xa2(k).ya2(k),xb2{k).yb2(k),Z2(k).R2a); 
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%  X-Y  Plot  of  the  estimate  of  the  coordinates  of  the  emitter 
figure(3) 


%  True  emitter  location 

subplot(1 1 1 ).  plot(xt/1  OOO.yt/1  OOO.**') 

hold  on 

%  Plots  the  final  locations  of  the  sensors 

plot(xa1  (k)/1  OOO.yal  (k)/1 000;r*'.xb1  (k)/1000.yb1  (k)/1  OOO/r*') 

plot(xa2(k)/1 000.ya2(k)/1 000.'r*'.xb2(k)/1 000.yb2(k)/1  OOO/r*’) 

%  A  priori  estimate  of  emitter  location 
plot(X0(1  )/1 000,X0(2)/1 000.'+’) 

%  Estimates  of  emitter  location 
plot(X(1  .:)/1000.X(2.:)/1000) 

%  Plot  three  sigma  error  elipsoids  for  various  points 
for  k  =  1:19:len-1, 

%  Selects  the  error  covariance  from  those  saved  in  vector  P. 
P1  =  [P(2*k-1,1)  P(2*k-1,2):  P(2*k,1)  P(2*k,2)]:' 

%  Calculates  the  coord,  of  ellipse. 
lxg.y1.y2,Emax,Emin]  =  ellip(P1,C,X(1,k),X{2,k)); 

%  Plots  on  graph 

plot((xg)/1 000.(y 1  )/1 000.'g-’,(xg)/1 000.(y2)/1 000, ’g-') 

%  Plots  the  loci  of  constant  TDOA 
plot(T1  (1  ,:)/1 000, T1  (2,:)/1 000,’:') 
plot(T2(1,:)/1000,T2(2,:)/1000,’:’) 

hold  off 

title(’Estimate  of  Position  of  the  Emitter:  tgt’) 

%  axis; 

axis([0  30  0  30]); 
pause; 


end 
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%  Plot  the  closeup  of  the  estimate  of  the  emitter  location 
figure(4) 

%  True  Emitter  location 

subplot(1 1 1 ).  plot(xt/1 000,yt/1  OOO,'*') 

hold  on 

%  Estirr  ?  of  emitter  location 
plot(X(  000.X(2,:)/1000); 

%  Calculate  and  plot  error  ellipsoids  for  steady  state  estimate 

[xg,y1,y2,Emax.Emin]  =  ellip(PK,C.xtavg.ytavg): 
plot((xg)/1000,(y1)/1000.'g-'.(xg)/1000.(y2)/1000,'g-‘) 

%  Plot  on  the  graph  the  major  and  minor  axis. 

Major  =  sprintf('%5.1f.Emax/1000); 

Minor  =  sprintf('%5.1f. Emin/1 000); 

text((xt-0.50*Emax)/1000,(yt+0.85*Emax)/1 000, [’Major  Axis  = '  Major '  km']) 
text{(xt-0.50*Emax)/1 000,(yt+0.75*Emax)/1 000, fMinor  Axis  = '  Minor '  km*]) 

%  Plot  the  loci  of  constant  TDOA 
plot(T1  (1  ,:)/1 000,T1  (2,:)/1 000,':') 
plot(T2(1,:)/1000,T2(2,;)/1000,’:') 

hold  off 

titleCCIoseup  of  Position  of  the  Emitter;’) 

%  axis; 

axis{[{xt-Emax)/1000  (xt+Emax)/1000  (yt-Emax)/1000  (yt+Emax)/T000]); 
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%  Ploci.m 
% 


% 

% 

% 

% 


% 

% 


% 


% 


This  program  calculates  the  coordinates  of  the  loci  of  constant  TDOA  for 
a  pair  of  sensors. 

Input:  The  coordinates  of  sensor  A  (meters):  xa,  ya 

The  coordinates  of  sensor  B  (meters):  xb,  yb 

The  Time  Difference  of  Arrival  (sec):  TDOA 

The  range  from  sensor  A  over  which 
the  loci  of  constant  TDOA  will  be  plotted:  R 

Output  The  vector  of  x  and  y  coordinates  of  the  loci 

of  constant  TDOA.  T 


% 

% 


function  T  =  ploci(xa,ya,xb.yb,TDOA,R); 


%  Calculate  the  distance  from  sensor  A  to  Sensor  B 
Rab  =  ones(1,length(R))*sqrt((xb-xa)''2  +  (yb-ya)'^2): 


%  Calculate  the  bearing  of  sensor  B  from  sensor  A 
psi  *  ones(1,length(R))*atan2(yb-ya,xb-xa); 


%  Calculate  the  difference  in  path  length  for  the  pulse 
D  =  3.0e+08*TDOA*ones(1.length(R)); 


Equation 


(|>  =  arccos 


[(R-D)2-R2-R^] 

-2R,bR 


+ 


V 


%  Calculate  the  bearing  required  for  each  range 

phi  =  acos(  ( (R-D).''2  -  (R).''2  -  Rab.''2  )./( -2*R.*Rab  ) )  +  psi; 


J  =  ones(1  ,length(R)); 

%  Calculate  the  x  and  y  coordinates 
T  =  [R.*cos(phi)+xa*J:  R.*sin(phi)+ya*J]; 
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%  ellip.m 
% 

%  This  program  calculates  the  coordinates  of  the  2D  error  ellipsoids  given  the 
%  error  covariance  matrix  and  the  location  of  the  estimate. 

% 


% 

input; 

The  2D  Error  Covariance  Matrix 

PK 

% 

The  number  of  standard  deviations 

C 

% 

% 

% 

The  X  and  y  coordinates  of  the  estimate 

xt.  yt 

Output 

The  X  vector  along  which  the  ellipsoid  is 

% 

% 

plotted. 

The  y  vectors  corresponding  to  the  upper 

xgout 

% 

and  lower  parts  of  the  ellipsoid 

y1out,  y2out 

% 

The  length  of  the  major  axis  of  the  ellipse 

Axmax 

% 

The  length  of  the  minor  axis  of  the  ellipse 

Axmin 

% 

function  [xgout,y1out,y2out,Axmax,Axmin]  =  ellip(PK.C,xt,yt) 


%  Inputs  the  error  covariance  matrix,  the  number  of 
%  standard  deviations,  and  the  location  where  plotted 

%  Calculate  the  inverse  of  the  error  covariance  matrix 
P  =  inv(PK): 

P11  =  P(1.1); 

P12  =  P(1.2): 

P21  =  P(2.1); 

P22  =  P(2,2): 

%  Calculate  the  upper  and  lower  limits  on  the  x  axis 

%  Equation 
f=sqrt(P11*C''2): 

%  Define  the  vector  along  the  x  axis  over  which  the  ellipsoid  lies 
xg  =  -f:(2*f)/1000:f; 
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%  Calculate  the  upper  ar^d  lower  curves  of  the  ellipse. 


%  Equation 


(PjlC^-JCtO 


y1  =  -(P12/P11)*xg  +  sqrt(  ((PirP22  -  P12''2)/Pir2)*... 
(P1 1*C*2*ones(1  ,length(xg))  -  xg.''2)  ); 

y2  =  -(P12/P11)*xg  -  sqrt(  ((P11*P22  -  P12''2)/P11''2)*... 

(P1 1*C*2*ones(1  ,length(xg))  -  xg.'^2)  ); 

%  Calculate  the  distance  of  the  ellipse  from  the  center 

Rellip  =5  sqrt(  xg.''2  +  y1  .''2 ); 

%  Calculate  the  major  and  minor  axis  of  the  ellipse 

Axmax  =  abs(2*max(Rellip)); 

Axmin  =  abs(2*min(Reliip)); 

%  Calculate  the  coordinates  of  the  ellipse  over  the  estimate 

y1out  =  y1+yt*ones(1,length(xg)); 

y2out  =  y2+yt*ones(1,length(xg)); 

xgout  =  xg+xt*ones(1,length(xg)); 

% 
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%  Burdist.m 
% 

%  Name;  Richard  W.  Wiliiamson 
%  Date:  18  February  1994 
%  Date  Revised;  18  February  1994 
% 

%  This  program  generates  a  simulated  burst  of  pulses  generated  by  a 
%  circularly  scanning  emitter.  The  length  of  the  burst  is  2.0  milliseconds  the 
%  maximum  amplitude  is  1.0,  and  the  PRF  varies  from  3  kHz  to  12  kHz. 

%  The  program  calculates  the  probability  density  function  for  the  centroid  of 
%  the  burst  for  a  peak  signal  to  noise  ratio  of  20  dB. 

% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

%  The  program  outputs  one  figure. 

% 

%  figure(1 )  A  plot  of  the  probability  density  function  for  the  burst  centroid. 

% 

%o^o^o^o/o^o^o^o^o^ozo^o^ozo^o^o^o^q/qzo^o^o^o^o^ozo^Q,^o^o^o^ozo^ozo^ozq/q/ 
/0/vA)/DA)/0/OAi/OAI/O/O/D/OAiA)/O/0/DAp/D/D/D/D/D  /D /D  /0/0/0f0A>/0/O/0f0/0 

% 

%  Input  Variables; 

% 

%  BW  =  Beam  width  in  degrees 
%  SR  =  Scan  rate  of  emitter  in  degrees/sec. 

%  PRF  =  Pulse  Repetition  Frequency  in  Pulses/sec 
%  PW  =  Pulse  width  in  seconds. 

%  FS  =  Sampling  Frequency  in  samples/sec 
% 

% 

%  Output  Variables; 

% 

%  Bur  =  The  noise  free  Burst  vector 
% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
clear; 

%  Initialize  the  variables 

BW=  1.00; 

SR  =  500.0; 

PRFk  =  [3000  6000  9000  12000] 

PW=  1.0e-06; 

FS  =  8.0e+06; 
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muz  =  {000  OJ; 
sigz  =  [000  0], 


zk  =  {0.90:(1.10-0.900)/200:1.10)/1000; 
dtz  =  zk(2)-zk(1); 

for  m  =  1:length(PRFk) 

PRF  =  PRFk(m); 

%  Calculated  Variables 


Bur_Len  =  BW/SR;  %  Burst  length  in  seconds; 

N_Bur  =  Bur_Len*FS;  %  Number  of  samples  per  burst 


N_PW  =  PWTS;  %  Number  of  samples  per  pulse  width 

N_PRF  =  (1/PRF)*FS:  %  Number  of  samples  per  PRF 


Bur  =  zeros(1  ,N_Bur);  %  Initialize  the  Burst  vector 
t  =  0:N_Bur-1 ;  %  Initialize  the  time  vector 


%  Load  the  individual  pulses  in  the  burst  vector. 

forl  =  1:N  PRF;N  Bur, 

Bur(l:l+N_PW-1)‘=  ones(1,N_PW); 
end; 


%  Define  burst  envelope  and  form  noise  free  burst  vector 

Bur_env  =  1 .0*sin(pi*l/N_Bur); 

Bur  =  Bur.*Bur_env; 


dt  =  1/FS; 

%  The  variance  of  the  noise  power  added  to  the  pulse 
V  =  0.010; 

ns  =  Bur; 

N  =  length(ns);  %  Length  of  the  burst. 
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%  Numerator  =  x  Denominator  =  y; 

%  Initialization  of  the  mean  and  sigma  variables 

mux  =  0; 

muy  =  0; 

sigx  =  0; 

sigy  =  0; 

cov  =  0; 


%  Calculate  the  statistics  of  the  numerator  and 
%  the  denominatior  of  the  centroid  equation. 


Np  =  0: 


fork=  1:N, 

Np  =  Np  +  1 : 
if  ns(k)  >  0; 

mux  =  mux  +  k*dt*ns(k); 
sigx  =  sigx  +  ((k*dt)''2)*V; 
muy  =  muy  +  ns(k); 
sigy  =  sigy  +  V; 
cov  =  cov  +  k*dt*V; 

end; 

end; 

sigx  =  sqrt(sigx); 
sigy  =  sqrt(sigy): 

rho  =  cov/(sigx*sigy)  %  The  correlation  coefficient 

c_bar  =  mux/muy 


dty  =  (10*sigy)/60; 

yk  =  (muy-5*sigy):dty:(muy+5*sigy); 

p  =  zeros(length(yk),length(2k)); 
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%  Calculate  the  probability  distribution  for  f(zy.y) 
for  I  =  1  :length(zk), 

for  k  =  1  ;length(yk), 

y  =  yk(k): 
z  =  zk(l): 

g1  =  (y*z-mux)/sigx; 
h1  =  (y-muy)/sigy; 

fxy  =  (1/(2*pi*sigx*sigy*sqrt(1-rho''2)))*... 

exp(  ( -1/(2*(1-rho''2))  )*(g1''2  -2*rho‘g1‘h1  +  h1''2)); 

p(k,l)  =  fxy; 

end; 

end; 


%  Initialize  the  variables, 
z  -  zeros(1,length(zk)); 

%  Integrate  the  probability  distribution  across  y  to 
%  calculate  the  marginal  density  for  z  the  centroid  variable. 

for  I  =  1:length(zk), 

H  =  0; 

for  k  =  1:length(yk)-1, 

H  =  H  +  dty*(  yk(k)*p(k.l)  +  yk(k+1)*p(k+1,l)  )/2; 

end; 
z(l)  =  H; 
end; 
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%  Check  the  summation  of  the  f(z)  density 
F  =  0: 

%  Change  the  scale  on  the  z  axis  scaled  In  terms  of  millisecond 
t  =  zk*1000; 


fork  =  1:length(zk)“1. 

F  =  F  +  dtz*(z(k)  +  z(k+1))/2: 

muz(m)  =  muz(m)  +  (t(k)+t(k+1  ))*dtz*(  z(k)  +  z(k+1)  )/4: 

end; 

F 

muz 

for  k  =  1  :length(t)-1 , 

sigz(m)  =  sigz(m)  +  ((t(k)-muz(m))'^2)*dtz*(  z(k)  +  z(k+1 )  )/2; 
end; 

sigz(m)  =  sqrt(sigz(m)) 

w  =  -0.5*(  ( t-muz(m)*ones(1,length(zk))  )/sigz(m)  ).''2; 

d_app  =  (1/(sqrt(2*pi)*sigz(m))rexp(w); 

dtt  =  t(2H(1); 

tgz(m.:)  =  z*dtz: 
tgapp(m,:)  =  d_app*dtt; 

end 

%  Output  the  probability  density  functions 

figure(1) 

plot(t,tgz/dtt) 

xIabeICCentroid  Time  of  Arrival  (milliseconds)’) 
ylabelfProbability  Density  (1 /milliseconds)') 
titleCSampling  Frequency  8  MHz') 
text(1.025,20,’3000  Hz') 
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text(1.015,32.'6000  Hz’) 
te)a(1.010.44;9000  Hz’) 
text(1.005.60.’12000  Hz’) 
text(0.96.60.’Peak  S/N  =  20  dB’) 
axis(I0.95  1.05  0  200.0]) 

%  sigzout  =  sprintf(’%5.4f,sigz(1)): 

%  muzout  =  sprinif(’%5.4f,muz(1)): 

%text(1.25,3,’S/N  15  dB’) 

%  text(1.25,2.5,rMean  = '  muzout]) 

%  text(1.25,2.0,[’Std  =  ’  sigzout]) 

%  text(-0.25, 5.0, 'Calculated  Density’) 

%  %  text(-0.25, 3.0, ’Gaussian  Approximation') 
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%  Puldist.m 
% 

%  Name:  Richard  W.  Williamson 
%  Date;  18  February  1994 
%  Date  Revised;  18  February  1994 
% 

%  This  program  generates  a  simulated  pulse.  The  length  of  the  pulse  is 
%  1 .0  microseconds  the  maximum  amplitude  is  1 .0,  and  the  sampling  rate  is 
%  8.0  MHz.  The  probability  density  function  is  calculated  for  peak  signal  to 
%  noise  ratios  of  15,  20,  25  and  30  dB. 

% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

%  The  program  outputs  one  figure. 

% 

%  figure(1 )  A  plot  of  the  probability  density  function  for  the  pulse  centroid. 
% 

%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% 

clear; 

%  Initialize  the  variables 
% 


muz  =  [000  0]; 
sigz  =  [000  0]; 

Vk  =  [0.03162  0.0100  0.003162  0.00100] 

%  initialize  the  variables. 

Dyk  =  [3;(9-3)/60;9;  3:(9-3)/60;9:  4;(8-4)/60:8;  5;(7-5)/60:7]; 

zk  =  3;(7-3)/200:7; 
z  =  zeros(length(Vk),length(zk)), 

for  m  =  1  ;length(Vk), 

yk  =  Dyk(m,;); 
dtz  =  zk(2)-zk(1); 
dty  =  yk(2)-yk(1); 


dt  =  1.00; 
% 
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%  The  variance  of  the  noise  power  added  to  the  pulse 
V  =  \/k(m); 

%  The  pulse  without  noise, 
s  =  [0  0.5  1  1  1  1  1  0.5  0 1; 

N  =  length(s);  %  Length  of  the  pulse. 

%  Generation  of  the  noise 
ns  =  s; 

%  Numerator  =  x  Denominator  =  y; 

%  Initialization  of  the  mean  and  sigma  variables 

mux  =  0; 

muy  =  0; 

sigx  =  0; 

sigy  =  0; 

cov  =  0; 

%  Calculate  the  statistics  of  the  numerator  and 
%  the  denominatior  of  the  centroid  equation. 

fork  =  1:N, 

mux  =  mux  +  k*dt*ns(k); 
sigx  s=  sigx  +  ((k*dt)''2)*V: 
muy  =  muy  +  ns(k); 
sigy  =  sigy  +  V; 
cov  =  cov  +  k*dt*\/; 


end; 

sigx  =  sqrt(sigx); 
sigy  =  sqrt(sigy): 

rho  =  cov/(sigx*sigy);  %  The  correlation  coefficient 
c_bar  =  mux/muy; 
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%  Calculate  the  probability  distribution  for  f(zy,y) 
for  I  =  1:length(zk), 

for  k=  1;length(yk), 

y  =  yk(k); 
z  =  zk(l): 

g1  =  (y*z-mux)/sigx: 
h1  =  (y-muy)/sigy; 

fxy  =  (1/(2*pi*sigx*sigy*sqrt(lH+io'^2)))*... 

exp(  ( -1/(2*(1-rho''2))  )*(g1''2  -2*rho*g1*h1  +  hr2)): 
p(k.l)  =  fxy; 

end; 

end; 

%  Integrate  the  probability  distribution  across  y  to 
%  calculate  the  marginal  density  for  z  the  centroid  variable. 

for  I  =  1  :length(zk), 

H=:0; 

for  k  =  1:length(yk)-1, 

H  =  H  +  dty*(yk(k)+yk(k+1))*(  p(k.l)  +  p(k+1.l))/4; 
end; 

2(m,l)  =  H; 
end; 

%  Check  the  summation  of  the  f(z)  density 
%  and  compute  the  mean  of  the  distribution 
F  =  0; 

%  Change  the  scale  on  the  z  axis 
t  =  (zk-ones(1,length(zk)))/(N-1); 

for  k  =  1:length(zk)-1, 

F  =  F  +  dtz*(  z(m,k)  +  z(m,k+1)  )/2; 

muz(m)  =  muz(m)  +  (t(k)+t(k+1))''dtz*(  z{m,k)  +  z(m,k+1)  )/4; 

end; 
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for  k  =  1;length(t)-1, 

sigz(m)  =  sigz(m)  +  ((t(k)-muz(m))'^2rclt2*(  2(m.k)  +  2(m,k+1)  )/2: 
end; 

sigz(m)  =  sqrt(sig2(m)) 
tg(m.:)  =  2(m,:): 


end; 

figure(1 ) 

t  =  (2k-ones{1,length(2k)))/(N-1); 
dtt  =  t(2)-t(1); 
plot(t,dt2*tg/dtt) 

title(’Sampling  Frequency  8  MHz') 
xlabel('Centroid  Time  of  Arrival  (microseconds)’) 
ylabel('Probability  Density  (1 /microseconds)') 

%  The  statistics  for  1 5  dB 
sigzout  =  sprintf('%5.5f,sigz(1)); 
muzout  =  sprintf('%5.4f,mu2(1)); 

text(0.545,25.’S/N  =  15  dB') 
text(0.545,20,['Mean  = '  muzout]) 
text(0.545,15,['Std  = '  sigzout]) 

%  The  statistics  for  20  dB 
sigzout  =  sprintf('%5.5f  ,sig2(2)); 
muzout  =  sprintf('%5.4f  ,mu2(2)); 

text(0.545,45,'S/N  =  20  dB') 
text(0.545,40,['Mean  = '  muzout]) 
text(0.545,35,['Std  =  ’  sigzout]) 

%  The  statistics  for  25  dB 
sigzout  =  sprintf('%5.5f  ,sigz(3)); 
muzout  =  sprintf(’%5.4f  ,muz(3)); 


text(0.545.65,’S/N  =  25  dB') 
text(0. 545,60,  [Mean  =  '  muzout]) 
text(0.545,55,['Std  = '  sigzout]) 
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%  The  statistics  for  30  dB 
sigzout  =  sprintf('%5.5f,sig2(4)); 
muzout  =  sprintf('%5.4f  ,muz(4)); 

text(0.545.85,'S/N  =  30  dB') 
text(0.545,80,['Mean  = '  muzout]) 
text(0.545,75,['Std  = '  sigzout]) 

%  axis([0.4  0.6  0  90]) 
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